{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "%matplotlib inline\n",
    "import sys\n",
    "sys.path.insert(0,'..') # allow us to format the book\n",
    "sys.path.insert(0,'../code') # allow us to format the book\n",
    "sys.path.insert(0,'./code') # allow us to format the book\n",
    "\n",
    "# use same formatting as rest of book so that the plots are\n",
    "# consistant with that look and feel.\n",
    "#import book_format\n",
    "#book_format.load_style('..')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook creates the animations for the Multivariate Kalman Filter chapter. It is not really intended to be a readable part of the book, but of course you are free to look at the source code, and even modify it. However, if you are interested in running your own animations, I'll point you to the examples subdirectory of the book, which contains a number of python scripts that you can run and modify from an IDE or the command line. This module saves the animations to GIF files, which is quite slow and not very interactive. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import filterpy.stats as stats\n",
    "import numpy as np\n",
    "from matplotlib.patches import Ellipse\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\lib\\site-packages\\matplotlib\\collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9e208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAEZCAYAAACad+BMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlwrNlZ3/Hvo33vTWurW9LdBrAheGwYTw0xnqIMmIkZ\nKolTMQFsTMV2nBhcUEWZIiY2IaRCAmUwIbaxjZfE2BATGzsMS6gwrmwMuDw23sZz79x7JXWrtfTe\n6n158ke3Bo2mJbV031a3up9PleqqpaN+z3u79dO7nHMeUVWMMeaogU53wBjTnSwcjDFNWTgYY5qy\ncDDGNGXhYIxpysLBGNOUhYMxpikLhx4hIu8UkXd0uh9nJSIfFpHXdbof5vksHHrHs6PZRGRERD4o\nIndFJC0iT4rIK0/6YRH5aRGJiEiq8bMjJ7T9JRH5soiUTwskEXGLyEdEZKfxcbS9Hu676R4WDr1p\nCNgAvltVZ4C3A78vIqvNGovI9wNvA74HWAWuAr94wvPfBH4W+CNO/8V+FzDWeN4HgB8TkR8/2oVT\nnsN0gIVDD1LVnKr+oqpuNB7/EXAHePExP/I64AOq+nVVTQL/GvjxE57/o6r6J0CG03+xXwX8B1Ut\nqOo68EHgJ860Q6YjLBz6gIgsAPcBXz2myQuALx16/DfAgoh4nOrCoc8HgG916HlNG1k49DgRGQY+\nBnxYVZ8+ptkUkDr0ON34d9qBLvwJ8DYRmRKR69SPGsYdeF7TZhYOPUxEBoD/DBSAt5zQdB+YOfTY\n1fg340A3fqqx/ZvAp4DfBcIOPK9pMwuHHiUiQv38fg74h6paPaH5V4EXHXr87cCOqiZa2NSJFyRV\nNaGqP6qqS6r6bcAg8EQLz2s6bKjTHTBt8x7gm4FXqGrxlLYfBT4sIh8DtoFfAD50XGMRGaL+3hkE\nhkVkDCipaq1J26vUT1mSwPcBbwC+++y7Yy6aHTn0oMYtyzdSPwLYFpFM4+OHG99faTwOAKjqnwL/\nHvgL4C7wDHDS+IUPADngNcC/bHz+o43nfpmIHD4deQn1C5xp4JeBf6KqX3dqX037iK0E1RsOBhep\n6knjE7qOiHwIeFxVP9LpvpjnsiMH0w3sL1QXsmsOveNxLucv2aepn8qYLmOnFcaYpuy0whjTlIVD\nj7Ap28ZpFg694znnhyLyFhH5vIgUGncEjiUioyLyLhEJi0hcRH6rMZahWds5Efl4o21SRP63iDxw\nyvP/iohEGx//rkm/7dy2C1k49K4w8EvA77TQ9ueoz9h8IfUJWi+mPs27mSnqIxxfDHiAjwB/JCKT\nzRqLyJuAHwL+TuPjBxtfe06zFvpoLpiFQ49S1U+p6h8CsRaavwr4TVVNqmoUeDfHTKtW1Tuq+uuq\nuqN17wdGqIdKM68DflVVt1R1C/hVTpgObrqHhUPva/Wv8tFp1QEROXVWpoi8iHo43DqmSbPp4C9s\nsU+mgywcel8r5/N/ArxVRGZFZJH6TEoFJk76IRGZoT7r852qetwMzmbTwada6JPpMBsE1ftaOXL4\nZcANfJH69OoPAC9S1Z1jn1RkHPgs8H9V9VdOeO5m08H3W+iT6TA7cuh9px45NJZw+0lVDajqdSAO\nfP649iIySn1k44aqHr24eFSz6eBfOb3bptPsyKFHicggMExjanXjF7rSbF0HEfE3Po0AL6V+p6Lp\nBcnGylKfpD4T88db6MpHgZ8RkceoH8X8DPAbZ9oZ0xGnHjmIyFsby5B/RUTeekybd4vITRH5kojc\n73w3zTn8AvVf4LdRn06dpz69+nlTtoFrwP+hfrj/IeBtqvrnxzzvQ8DfA74XSB6aDv5djed+zpRt\nVX0f9dOPL1O/GPlZVf1tZ3fVtMOJcytE5FuBjwPfCZSpX7j6Z6r6zKE2jwBvUdVHROSlwG+o6oPt\n7bY5yqZsG6edduTwzcATjXPSKvA54B8cafMo9YEwqOoTgLux2rExrbIRkl3otHD4CvAyEfGKyAT1\nw8nAkTbLwOahx6EmbUz7PU59JafL5tM8dxyE6RInXpBU1adE5FeAPwOywJPA89YJ5Pm3y+wvwQVT\n1c91ug/n0RjFabrQqXcrVPV3aIzPF5F/S73M2mFhIHjocYAmS4+LiAWGMR2iqmeev9LK3Yr5xr8r\nwN+nXnfgsM8Ar220eRBIHjd4RlUv7OMd73iHbe8Sbsu25/zHebUyzuGTIuKjfrfin6tq+mBWnaq+\nT1UfE5FHROQW9VOP15+7N8aYrtHKacXzagxo/d714ccnVVMyxlxCPTt8+uGHH7btXcJt2fa6x4Ut\nMCsielHbMsb8LRFB23FB0hjTnywcjDFNWTgYY5qycDDGNGXhYIxpysLBGNOUhYMxpikLB2NMUxYO\nxpimLByMMU1ZOBhjmrJwMMY0ZeFgjGnKwsEY05SFgzGmKQsHY0xTFg7GmKZaWX36pxt1Mr8sIr/b\nKMh6+PsPi0hKRJ5sfLy9fd01xlyUExeYFZFl4CeBb1HVooj8HvAaGuXvDvmcqj7apj4aYzqglaXp\nh4AJEakCEzQpWMPzK14ZYy65E08rVDUM/Br1Kldb1AvWHC3NrsBDIvIlEXlMRF7Qnq4aYy7SieEg\nIh7qVbTXAD8wJSI/cqTZF4Cgqn478JvUC6MaYy65004rXgHcUdUYgIj8N+Ah4GMHDVQ1c+jzPxaR\n/yQiXlWNH32yd77znc9+/vDDD1+a9fuNuUwef/xxHn/88Xt+nhPrVojIA9SL6H4nUAA+DPyVqv7W\noTYLwK6qaqP976vqWpPnsroVxnTAeetWnHjkoKp/JSKfpH7qUGn8+/7DtTKBVwNvFpEKkKN+N8MY\nc8lZxStjHFAqV9nPl8gWShTLVUaGBpkcG8bnmuh019pz5GCMOZ6qktwvsJfMkcjmyJVy5Ct5ytUS\nQwPDuMdcrBR9BOddne7quVg4GHNGqspeMsd2fJ9kPkOikCBXzjIxNsDY2CATw4MUyzk2UimGBodZ\n8EwxMjzY6W6fmYWDMWcQT+fZimWI7SfZze0hA2Xc0yMsTU0xOPC3R+7TDLOfq1CqliiWKxYOxvSq\ndLZIOJomtp9mN7tHlQLzvjGmJsaO/ZmBAajVqlRrl/Nam4WDMScoliqs76SIZjLsZXcpVLPMesZw\nT0+f+rOFYo2xiXEmRocvoKfOs3Aw5hi7iSybe0m2M7tkykm8rlGWZ6YROf3Cf7FUZYBBxodHLuUp\nBVg4GPM8pXKVu9tJdtJJtjJbTIzD1YXp51xTOE0yU2Jm1IVrcvT0xl3KwsGYQ/aSWUJ7aSKZHTLl\nJIuzY0xNnO20QFVJZ8uszriZ7YJxDudl4WAM9aOF9Z0kO6kkW5kIY2PKmn+SocGzL5aWzpYZHZjA\nPTnB5PhIG3p7MSwcTN+LpnJs7qbYzuySKsVZ9I0zPXm+i4iqyl6iwNJkkHn3pMM9vVgWDqZvlStV\n1ndS7CSTbO1HGBmtcmV56lxHCweSmRKjAxP4Jme6Yuj0vbBwMH0pkcnXLzru75IsxlnwjTEzefyY\nhVbUakosVSIwtcry3IxDPe0cCwfTV1SV0F6aUDRJOBNieKTKFf8UQ0P3vhB7PFVkfHCK2elp3FP3\nFjTdwMLB9I1SucrtSILtVIzIfoQ5zwieGWeuC1SqNRLpMmvuQE8cNYCFg+kTyf0CdyIJttLb7FeS\nrCxOMjbq3OCkaKLIzKib2Zlppi7xHYrDLBxMT1NVwtEMoWiCcDrM0EiFK8tnG9B0mv1cmf1sjWve\nOQI9ctQAFg6mh5UrVW5vJdhOxYnsR/C5h/G6nL29WK0p29E8S9NBAnMuxkZ651eqd/bEmEP28yVu\nb9VPI1KlBIGFccbHnH+778byTA27mJ9xs+idcvz5O8nCwfSc3USW9Z0E4UyY2kDh3CMdT7OfK5PN\nw3XvAmuLbsefv9PuuVZmo827ReRmo7DN/e3pqjEnq9WUO5EEN7d2uJ28w+h4mdWlexvUdJxqTdmO\nFViaWiIw21unEwdOK2pzUCvzJar6bcAgR1aXFpFHgOuqegN4I/CeNvXVmGMVSxWe2ohye3ebzfQG\nC75h5r3jbdveTizP1JCLhRk3Cz12OnHAiVqZj9IorKuqT4iIW0QWVHXH2a4a01wmV+SZrTjhVIRs\nNc3K0gSjI+1bQ2E/VyaXh+ve+Z48nTjgRK3MZWDz0OMQEHCyk8YcJ5rK8dTGHrfj61QG9rnin2pr\nMBycTvin/QTn3Iz24OnEgRP37EitzBTwX0XkR1T1Y0ebHnncdNE8K4dnnBTaS7OxF2czvcnMlDDv\nbf8syO1onulhNwszbuY93Tnr8qLK4f0j4PtV9Z82Hv8Y8KCq/otDbd4LPK6qn2g8fgp4+dHTCitq\nY5xSrda4s51kKxFnKxNmzjuCe7r9oxLT2RJ7sSrXvVd54dr8pTlqOG9Rm9Mu464DD4rIuNQXznsF\n8LUjbT4DvLbRiQepn3rY9QbTFqVylW9sxrgbjbC1v8nywtiFBEOpXGUnWmR52k9w3nVpguFe3HOt\nTFV9TEQeEZFbQBZ4fbs7bfpTNl/iVjhOKBUhV02x5p9i2IHZlKdRVcK7OWYn5ljyeJi75Iu4tMpq\nZZpLIZ7Oc2c7zmYqBINF/PMTjs6POMl2NE+1NMZV3wrfsjLLYBvGTbST1co0PWsrmmlceAwxMaEs\n+i7uL3dqv0Q2B9e8S1zzey5dMNwLCwfTtWo15e52knAiTjgdYtYzjGfm4hZRKZaq7MaKBGdWWV3w\nMH5Ji9Ocl4WD6UrlSpVb4ThbySh7uR388+NMjl/c27VWq19nmJuYx++93EvMn5eFg+k6uUKZW+E4\n4dQ26XKclaXJtg5saiYSzTM+OM3izCyrC707CvIkFg6mqyT3CzyzFSOUClMdyDu+MEsr4qkipeIA\n17x+ri17Gbjg7XcLCwfTNfaSWe5sx9lIbTA2VmN5drKlupROyhUqxJJl1txrrC26e3K2Zav6d89N\nV9mKZljfjbOR3sA9M8Cs++LP8SuVGlt7eZamlwnMevBMt29W52Vg4WA6bn07SShWnyNxUUOhj1JV\nQrs53CNellzenloL8rwsHEzH1GrK7UiCcCJGJBNmae7sRWudshMvMKTj+F0LXPV7OtKHbmPhYDqi\nUq3V70gk9tjLbxNYmGjLGo+tSGZKZLNwxe3n6pKnLStHXUYWDubClcpVboZihFI7JItRVpcmGRm+\n2FuVB3KFCnvxEquuNa4seS91VWynWTiYC5UvlrkZirOZ3CJfS7G25EwpuvMolauEd3P4p4ME5zx9\nOdDpJBYO5sJkckVuheNspELUBnKsLE1d+BiGA9WasrmTY258gWWPXYBsxsLBXIhEJs+tcIzNdIih\nkTLBuYsfw3CgPgU7y9SQmyXXHFeW7AJkMxYOpu12E1nubMfYSG8yOaEs+Dp7+L4TKyDVMQLeJa73\n8QjI01g4mLYK76VZ34uzmdrA4x7C5+psafpEukguL1z1BLjm93TsQuhlYOFg2kJVWd9JEYrFCKVD\nzPtGcE119k7Afq5MLFlh1bXKlUW7M3EaCwfjuFpNeWYrzlYiRmR/q6ODmw4US1UiewUCM0FW5rx4\nZ/p7aHQrLByMo6rVGjfDccKJXaKFXYKLE4yNdvbQvVKtsbmTZX5yCb/Hi392uqP9uSxaqZX5TSLy\n5KGPlIj81JE2Dze+ftDm7e3rsulWlWqNp0MxNuLbRAs7rHRBMKgqoZ0crhEvfvdsT1eoctqpRw6q\n+g3gfgARGaBeDu9TTZp+TlUfdbZ75rIoV6o8vRkjlNwmVY53dHDTYZFoniEmWHYtcs3vsTsTZ3DW\n04pXAM+o6maT79n/ep8qlircDMfZSGyRqyZZXWpPyfuziiYLlAqDrHn8XF/2MjxkdybO4qyv4GuA\n323ydQUeEpEvichjIvKCe++auQwKpQrf2IyxHg+Rr6VYaVPJ+7PKZMskU1UCMwGu+X19tzisE1o+\nchCREeAHgbc1+fYXgKCq5kTkB4BPA/cdbWS1MntLrlDmZijGRipMmX1WFie74rA9V6iwHS0QnFlh\ndcGLe6qzYysu2oXUynxOQ5EfAt6sqq9soe0d4CWqGj/0NStq00Oy+VL94mNyEx0sEFiY6Nhw6MOK\npSob21n8U0HW5uZYWXB1uksd165amYf9MPDxYza+0KiliYg8QD104s3amssvkyvy1OYedxMbMNw9\nwVCp1NjcyTE/scSyx2vBcI9aOq0QkUnqFyPfcOhrz9bLBF4NvFlEKkCO+rUJ04NS+wVubcW4m9hg\ndKzC0lx31I2s1pSN7Sye0VmW3bO2mpMDrFamadnBzMqN1AZjEzUWfd0xylC1Hgxj4mLNu8w3BX19\nVbbuNFYr07RVNJXjdiTGemqD6SmY83RHMACE93IM6SRBzxI3Al4LBodYOJhT7Say3N6Osp46WDa+\ne67+b8fy1EqjrHmWuRHw2VgGB1k4mBNFYhnu7sTZSK3j8wzhmRntdJeeFU0WyOeEK54ANwK+vi5A\n0w72v2mOFd5Lc3c3Riizyax3uCP1JI6TzJRIpmtcca9xze9jyqZfO87CwTS1sZNiIxollA6xMDvC\nzGT3/PLt58pEE2VWZurrMvTbIKeLYuFgnmd9O8lGLEY4vdkVazEcli9UiOw1Rj/Oe5lzd8et1F5k\n4WCe43Aw+OfHmRzvnrdIqVwltJtjaTpAwGfrMrSb3fMxz+rmYKhUamxsZ5kbX2TZ42PV1mVou+55\n9U1H3d1OshGNEs6ECCyMM9Gh0nTN1GtMZHGPzrLsmeOqLSV/IezIwXR1MNRXcsoyPuAi4FqwpeQv\nkIVDnzsIhq0uDYbwbmP0o9vPjYCvK9aK6Bfd804wF+5OJMFGLEokE2a5y4IBYGsvD5UxVj0B7gv6\nrMbEBbMY7lOHg6HbjhgAIns5qqURVtxBG/3YIRYOfehoMIx3WTDsxPIUC0OsuOvDoifGumecRT+x\ncOgzd7eTbMZiXRsMe4kCuZyw6l7hxvKsDYvuIAuHPnL04mO3BUMsVSSzDyuuFa4v+5iZ7J5JXv3I\nwqFPHA6G5S4MhkS6SCJVZdW1wnX/rM2X6AIWDn1g/dA4hm68K5HaLxFLVFh1rXB1yWd1LLuEhUOP\ne3ZIdBeOY4B6fYm9eJmga4UrCz5mXROd7pJpODEcWqmT2Wj3bhG52Shqc3/7umvO4vBcieX57guG\n/VyZ7WiRwHSQKws+FrxTne6SOeTEd0srdTJF5BHguqreEJGXAu8BHmxPd02rNndTXTuJCuqFZyJ7\nBQIzQVbnfSz5bIZltznLacVxdTIfBT4CoKpPAG4RWXCof+YcwntpNvacC4b3ffbzDvWsrlCsEt7J\n458OsDI7S2BuxtHnN844SzgcVydzGTgcGCEgcC+dMue3Hd9nfS/OZmOhFieOGN7/37/gQM/qCqUq\nmztZlqb8BH0+KzzTxVotanNSnUx4foXtpgUqrFZme+0mstzZjrGZ3mBhduSeV3B632c//2wwfMeb\nfps3vOrFvOkHv+Pcz1cqVwnt5FiY8BPwzrJmazK0xYXWyjypTqaIvBd4XFU/0Xj8FPByVd050s6K\n2rRRLJXj5tYe68l1xxeD/Y43/Taff98b7+k5ypUa65Ess2MLrPjmub7s7YoSev2g3bUyj62TCXwG\neG2jEw8CyaPBYNorkcnzTCTGRmoTr3vI8VWi3/CqF9/Tz1cqNdYj+/jG5gh45rjmt2C4DE49cmjU\nyVwHrqhqpvG1w3UyEZH/CLwSyAKvV9XnnaTakUN7pPYLPB2Kcje1zswUzHq6a2RhpVLjbmQf79g8\nQc8CN5atItVFO++Rg9XKvMQyuWI9GBLrTEzWmPd218hCC4bu0O7TCtNlsvkSt8Jx1pObjI13ZzCs\nb2fxjM4RcM9bMFxC9mpdQvlimadDMdaTGwyPllmc7c5gcI/MEvQscF/Aql5fRvaKXTKFUoWnN2Ns\npELIcBH/XHfNRahU68HgGvFZMFxy9qpdIqVylZuhGBupMDqQZ7kbgyFSD4YVzyI3LBguNXvlLoly\npcrToRgbiTBl9lmen+iq24GVao2NSBbXsJege9FWiu4B9updApVqjac3Y2wmIhQ0Q3BhsqtqNxwE\nw8ywl6BnifuCFgy9wF7BLlet1rgZihFKbrNfSRBc7LJgaIx8tGDoPd01j9c8R62m3ArHCSV3SZXj\nrC5NMthtwXBw8dG9aMHQY+yV7FKqyu1IglBij3h+l5XFya76xTscDCseC4ZeZEcOXWpjJ0U4EWMv\nv83q0iTDQ93zi2fjGPqDvaJdaCuaIRSLs5UJE5if6KoycAdDoi0Yep8dOXSZvWSW9d36Yi3+ubGu\nWkK+3Jhd6R2bJ+Cet2Docd3zzjMk9wvc2Y6zkd5gznvvi7U46XAw2CSq/mDh0CX28yWe2YqxkdrA\nPTPg+JoM96JYqrK5k8M3Nk/AgqFvWDh0gXyxXB8WnQwxPqHMurtnIlWhVF/abW58kYBnjuvL3q4a\nZ2Hax8Khw0rlKrfCcTZTYQZGiiz6JjvdpWflCxXCu3kWJpdY9sxyzW/B0E8sHDqoWq3VgyEZoSJZ\nVua6Jxiy+Qpbu3mWppcJeme5suTuqrkcpv0sHDpE9WD04w7ZSpJV/1TX/PLt58pE9goszwQJ+ny2\nSnSfsnDokNtbCcLJKIlilLUuGhadzpbYiZYIzqwQnLW6Ev3s1EvOIuIWkU+KyNdF5GuNFaYPf//h\nRg3Ng3qab29fd3tDffRjnN1chODCBENdMvoxmSmxGysTnFlhbX7WgqHPtXLk8BvAY6r6ahEZApqd\nGH9OVR91tmu9KRLLsBmNsZUJEViYYHSkO0Y/xlNF4skKK65V1ha8VrvSnBwOIuICXqaqrwNQ1QqQ\nata0DX3rOdFUjvXdOKFMiMXZ7hn9GE0WSKZrrLpXubo4y7yney6Mms457Xj2CrAnIh8SkS+IyPtF\n5OjaZAo8JCJfEpHHROQF7enq5ZbaL3C7UXhm1jPM9GR3jH7cjedJZ+CKe43r/jkLBvOs0/50DQEv\nBt6iqn8tIr8O/Bzwrw61+QIQVNWciPwA8GngvmZP1q+1MvPFMrcjCTZTm8xMS9eMftyO5SnkBlhz\nrXB92YdnunsGX5nzu5BamSKyCPw/Vb3SePx3gZ9T1Ved8DN3gJeoavzI1/uyqE2lWuPr63vciW/C\nUJ7l+e5YFDayl6NYGGLNs8J1vw/XVHdVyjLOaUtRG1XdBjZF5OBI4BXAV49seEEaN+hF5AHqgRPH\noKo8E46zldqhTBb/XOf/Mqsqod0s5eIIa54V7gvMWjCYplq5IvaTwMdEZAR4BviJI7UyXw28WUQq\nQA54Tbs6e9nc3U6ylYyRLMVZ8092fJCTqhLaySHVMdY8QW4EfEyNd8cpjuk+ViuzTSKxDLe399hM\nbxBcmmCsw7csqzUltJNlSCdZdQe4EfAxMdYdF0VNe533tKI77qX1mEQmz8ZuglAmxNLcWMeDoVKp\nsbGdZXLITdCzxI2Aj7ERe+nNyewd4rBcocztSH0lJ49rqOMLthysxeAZnSXgXuD6srerlp0z3cvC\nwUHlSn36dSgVZmyshs/V2TsTuUKF8G6O+Ykllt2zXLdFWswZWDg45KDGRDi1TXUgz/JsZwcTZbJl\ntqMF/NNBlj1ervo9Hb8gai4XCweH3Ikk2EpGyVQSrHV4+nUyUyKaKNvMSnNPLBwcEN5Ls5VMPFtj\nopPTr/cSBdIZZdW1yuq8TaAy52fhcI9iqRwbewnC6RD+DteYqI96HGTNvcLVJR+zHb7mYS43C4d7\nsJ8vcTtSn2U55x1hcrwz/521mrK1l0Mro6x5AjYc2jjCwuGcSuUqz4TjhNJhJia0Y5OpqjVlczvL\niEwR9Pi5sexj0kY9GgdYOJzDQZHbcHobGercitHlxuCm6SEPAfeiDW4yjrJ30jls7KSIJGNkygmu\nLHfmgt9BPQnv2BzLrnluBLwMD9ngJuMcC4cziqZyhONJdrIRgh26M5ErVAjvHK4n4bHBTcZxFg5n\nkCuUWd9JEEqHmPOOdmTORH116CL+6QABr8/qSZi2sXBoUaVa45mtOKH0FpMTdOQCZDRRIJmpEZxZ\nZWXOS3DeBjeZ9rFwaNGdSIJIeo8KOZYv+AKkqrK1l6dcHOKKe4W1Ba+t9WjazsKhBZFYhkgy8WwB\nmos8jK9UaoR2c4zIFFe9fq75vcxMjl7Y9k3/snA4RSZXZHMvSWR/i6W58QstQHNwR8I96mPZVZ9u\nbbcqzUWxd9oJypUqdyJJQukw7unBCx0BmcmWiUQLLE0tseTycW3Zy5DdkTAXyMLhBHciSbYaA51m\nPVMXtt1YqkgiVWFlZoVln5fVBZfdkTAX7p5rZTbavFtEbjYK29zfnq5erK1ohkgqTqqUwD93MROY\nVJXIXo50WllzrXFtaZ61RbtVaTrjnmtlisgjwHVVvSEiLwXeAzwvQC6TdLbI5l6C7f0I/vnxCzmc\nr1RrhHdzDOoEVzx+rvl9uG3ylOkgJ2plPgp8pPH9JxpHGguqutOODrdbtVrj7naScGYL98wgExdQ\nz7K+zmMW14iXZdci1/xeWxnadJwTtTKXgc1Dj0NAwMlOXqSN3RQ7+1EYLDLrbv9f7v1cmY1Ijrnx\nJdZ8y3zzyqwFg+kKTtTKhOdX2W5aoKLba2Um9wtsJ9JEc3us+dt/nSGRLhJLVgjMBPF7vKwtuhno\n4CpSpjd0Ta1MEXkv8LiqfqLx+Cng5UdPK7q9qE2lWuNrd/e4Gb2Ny6V4Zto30EhV2Y7myRcGCM4E\nWJnz4p+15dxMe3SsVibwGeC1jU48CCQv4/WGjZ0U25k9ZKjU1mAoV2qsR7LUyuNc9axx3/K8BYPp\nSvdcK1NVHxORR0TkFpAFXt++7rZHIpNnO5kino+ytty+OQv7uTKRaB7f+DxLM3Nc83sYH7XrC6Y7\n9X2tzEq1xlfv7HIrdhu3W9o22zKaLJBMV/FPL7Pk9nBl0W1rMJgLYbUyz2l9O8n2/h4Dw2Xc086P\ngqzW6gObKqVhrrhXCM65bbl4cyn0dTjE0/XTiUQ+xpU2nE4UilVCu1mmhtys+ha5umQzKs3l0bfh\nUK5U2dgI69jKAAAISUlEQVStz7ac8446PtsymSmxFy+yMOln0eXl6pKHUZtRaS6Rvn23hqMZdvaj\nDA1XcU87N9ipVlMi0Tyl4gCrriss+9wE52dsfoS5dPoyHPLFMjvJDIl8zNG7E4VSlfBujsnBGa55\nl1hbdOOZHnfs+Y25SH0ZDuFohmg2yszUEMMOnU4cnEbMTy6yOOPjqt9jC7OYS63v3r2ZXJG9VIZ0\nKcnV+Xu/O1Gp1NiK5qiVh1lxrbHsdROcd9kwaHPp9V04hKMZ9nJ7eF0j9zwV++BowTPuY9E3x8qC\ny04jTM/oq3BI7heIZdLkq/sszZz/qKFSqRGJ5qmUBwnOrLLodrG66LZl3ExP6atwiMTqRw2z7tFz\nH/an9kvsxgt4Rn0s+OZYmXfjnbGjBdN7+iYccoUyqVyOQjVHYOrsIxQr1Rrb0Tyl0iDB6TUW3C5W\nF1xWn9L0rL4Jh71klkQ+yczk8JnHHKSzJXZiBdyjXla886zMu/C5LmZdSWM6pS/CQVVJ7BdIFZOs\neFs/BSiVq+zGC42jhVXmXS5WF9yMDNvRgul9fREOhVKFXKmADNQYbaH4baVSI5oskslW8Y77WPH6\nCM67mLWjBdNH+iIcsoUy+XKe8dGTd7daU+KpIol0Gfeoh2seHwueaZZ8U3ZtwfSdvgiHarVGVasM\nDTa/1lCp1Ejul0ikS0wNz3DVHWTONc3y7LRNljJ9qy/e+SPDg4wMDpMpVqnVlIEBoVKpkS9WSWVL\n5PI1ZkZmWJlZZna6HgqT4+1Z9MWYy6IvwmFqfAT3+AypYpqn19P1MQ46wPjwGNMjsyx7XXinJ5hz\nTzA9YestGAN9tExcrlDm7naSQqlCuVZhdGiYybFhpidG8c2M2zUF07POu0xcS+EgIneBNFAFyqr6\nwJHvPwz8IXC78aU/UNV/c6RNV6whqapUqjULA9M32r2GpAIPq2r8hDafU9VHz9qBiyYiFgzGtOAs\nM4VOSx6bo2xMD2k1HBT4cxH5vIi84ZjvPyQiXxKRx0TkBc510RjTCa2eVnyXqkZEZA74HyLylKr+\nr0Pf/wIQVNWciPwA8GngvqNP0u21Mo3pBRdSK7PpD4i8A9hX1V87oc0d4CWHr1F0ywVJY/pNW2pl\nNp54QkSmG59PAt8HfPlImwVpTHUUkQeoh85JFy+NMV2uldOKBeBTjd/9IeBjqvpnh+tlAq8G3iwi\nFSAHvKZN/TXGXJC+GQRlTL9q22mFMaY/WTgYY5qycDDGNGXhYIxpysLBGNOUhYMxpikLB2NMUxYO\nxpimLByMMU1ZOBhjmrJwMMY0ZeFgjGnKwsEY05SFgzGmKQsHY0xTFg7GmKYsHIwxTVk4GGOaamWB\n2bsi8jci8qSI/NUxbd4tIjcbdSvud76bxpiL1sqRw0EpvPuP1sgEEJFHgOuqegN4I/Aeh/t4Lk6s\n22/bu/ht2fa6R6unFSctTvko8BEAVX0CcIvIwr127F71+gtu4WDba7dWjxxOKoW3DGweehwCAk50\nzhjTOa3UrTitFB48/8jC1qA35pI7U92KZqXwROS9wOOq+onG46eAl6vqzpGftcAwpkPOU7fixCMH\nEZkABlU1c6gU3i8eafYZ4C3AJ0TkQSB5NBjO2zljTOecdlpxaik8VX1MRB4RkVtAFnh9W3tsjLkQ\nF1YOzxhzuTg+QlJEXikiTzUGRb3tmDaODJo6bVsi8rCIpBoDuJ4Ukbefd1uN5/sdEdkRkS+f0Max\nAWGnbc/J/RORoIj8hYh8VUS+IiI/dUw7p167U7fn8P6NicgTIvLFxvbeeUw7p/bv1O05/f5sPOdg\n47k+e8z3W98/VXXsAxgEbgFrwDDwReBbjrR5BHis8flLgb9s47YeBj7j4P69DLgf+PIx33dk386w\nPcf2D1gEXtT4fAr4RrteuzNsz+nXb6Lx7xDwl8BL2/z6nbY9R/ev8Zw/A3ys2fOedf+cPnJ4ALil\nqndVtQx8AvihI22cGjTVyrbg5AFcZ6L1W7iJE5o4OiCshe2BQ/unqtuq+sXG5/vA1wH/kWaO7V+L\n2wNnX79c49MR6n9QakeaOP36nbY9cHD/RCRAPQA+cMzznmn/nA6HZgOilltoc55BU61sS4GHGodQ\nj4nIC86xnXvtUzsHhLVl/0RkjfoRyxNHvtWW/Tthe47un4gMiMgXgR3gz1T1r480cXT/Wtie06/f\nu4CfpXkIwRn3z+lwaPXqphODplr5mS8AQVX9duA3gU+fYztndZEDwhzfPxGZAj4JvLXxF/15TY48\nvqf9O2V7ju6fqtZU9UXUfyFeKiIvbNaloz/Wxu05tn8i8ipgV1Wf5OSjkZb3z+lwCAPBQ4+D1NPp\npDaBxtcc35aqZg4O7VT1j4FhEfGeY1vn7dN5960lTu+fiAwDfwD8F1Vt9kZ1dP9O2167Xj9VTQF/\nAbzyyLfa8vodtz2H9+8h4FERuQN8HPgeEfnokTZn2j+nw+HzwA0RWROREeAfUx8kddhngNcCnDRo\nyolticiCNAZpiMgD1G/dxs+xrVY5tW8tcXL/Gs/zQeBrqvrrxzRzbP9a2Z7D+zcrIu7G5+PA91K/\nznGYk/t36vac3D9V/XlVDarqFeA1wP9U1dceaXam/WtlbsVZOlgRkbcAf0r9bsIHVfXr0oZBU61s\nC3g18GYRqQA56v9p5yYiHwdeDsyKyCbwDuoXmhzdt1a3h7P7913AjwJ/IyJPNr7288DKwfYc3r9T\nt4ez+7cEfEREBqn/Ufy9xv60a0DfqdvD4ffnEQpwL/tng6CMMU3ZMnHGmKYsHIwxTVk4GGOasnAw\nxjRl4WCMacrCwRjTlIWDMaYpCwdjTFP/H99lRy0GALXOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9e320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "#import mkf_internal\n",
    "import matplotlib.pyplot as plt\n",
    "from gif_animate import animate\n",
    "\n",
    "def ellipse_animate(frame):\n",
    "    \n",
    "    plt.cla()\n",
    "    cov = frame/15.\n",
    "    P = np.array([[2,cov],[cov,2]])\n",
    "    stats.plot_covariance_ellipse((2,7), cov=P, facecolor='g', alpha=0.2, \n",
    "                              title='|2.0 {:.1f}|\\n|{:.1f} 2.0|'.format(cov, cov))\n",
    "fig = plt.figure()\n",
    "animate('multivariate_ellipse.gif', ellipse_animate, 30, 125, figsize=(4, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='multivariate_ellipse.gif'>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.0025,  0.005 ],\n",
       "       [ 0.005 ,  0.01  ]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f.Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\lib\\site-packages\\matplotlib\\collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFpCAYAAACReS/IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VNX5wPHvnSXrTCb7zpKw71tQQAVEUVGLUq1t1SrW\nirVqtdYNRYOKKJaibaXan9YNqyJarVqKVBGUshgQBMKOAUI2siczyWS2+/tjSDKTfc+QeT/Pk8d7\nZ849970mzDv33LMoqqqqCCGE8Fua3g5ACCFE75JEIIQQfk4SgRBC+DlJBEII4eckEQghhJ/T9dSJ\nysvLe+pUQgjht0wmU7uPkTsCIYTwc5IIhBDCz/VY05Cnjty6dJcdO3YAkJaW1suRNOarsUlc7eOr\ncYHvxiZxtU9nm97ljkAIIfycJAIhhPBzkgiEEMLPSSIQQgg/J4lACCH8nCQCIYTwc5IIhBDCz0ki\nEEIIPyeJQAgh/JwkAiGE8HOSCIQQws9JIhBCCD8niUAIIfycJAIhhPBzkgiEEMLPSSIQQgg/J4lA\nCCH8nCQCIYTwc5IIhBDCz0kiEEIIPyeJQAgh/JwkAiGE8HMtJoKVK1cybtw4TCYTJpOJadOmsXbt\n2rr358+fj0aj8fqZNm1atwcthBCi6+haerNfv34899xzDBkyBJfLxRtvvMHVV1/Nzp07GTNmDIqi\nMHv2bFatWlV3TEBAQLcHLYQQouu0mAjmzp3rtb9kyRJeeukltm3bxpgxY1BVlYCAAGJjY7s1SCGE\nEN2nzc8InE4n7733HhaLpa75R1EUNm/eTFxcHMOGDWPBggUUFhZ2W7BCCCG6nqKqqtpSgb179zJ1\n6lRqamowGAy88847zJkzB4DVq1cTGhpKSkoKWVlZLFq0CKfTyc6dOxs1EZWXl9dtHzlypBsuRQjR\nF0SuXUtAYSE1yclUjh+PIyqqt0PyeUOGDKnbNplM7T6+1URgt9vJzs6mvLycNWvW8Morr7Bx40ZG\njRrVqGxeXh4DBgxg9erVzJs3z+s9SQRCiLYY+pvfEJaRAcCRFSsov+CCXo7I93U2EbT4jABAr9eT\nmpoKwIQJE8jIyOD555/n1VdfbVQ2ISGB5ORkjh492mKdaWlp7Q60u+zYsQPwrZhq+WpsElf7+Gpc\n4KOxnT5dtzlkzhwYPrwXg/Hmk/+/8P6i3RGtJoKGnE4nNputyfcKCwvJyckhISGhU0EJIfyUqsKP\nfkT5jh0EnjpFUEpKb0fkF1pMBA8//DBXXnklycnJVFZW8s4777Bp0ybWrl2LxWIhPT2da6+9lvj4\neI4fP87ChQuJi4tr1CwkhBBtoijwl79wpPabd2BgLwfkH1pMBAUFBdx4443k5+djMpkYN24c69at\nY/bs2VitVvbt28eqVasoKysjISGBWbNm8cEHHxAaGtpT8QshhOikFhPB66+/3ux7QUFBrFu3rssD\nEkII0bNkriEhhPBzkgiEEMLPtbvXkBBCdJsXXoDoaELtdiwjRvR2NH5DEoEQwjfY7fDAA+BwMAL4\n7ptvejsivyFNQ0II33DyJDgcANhiY3EFBfVyQP5DEoEQwjd4zEhQk5zci4H4H0kEQgjf4JEIrP36\n9WIg/kcSgRDCN3jeESQl9WIg/kceFgshfMPFF7v/e/QoVT400Zw/kEQghPANV1zh/gEqzsw1JHqG\nNA0JIYSfk0QghBB+ThKBEEL4OUkEQgjh5+RhsRCi9/3jH5CZCYMHw6xZvR2N35FEIITofR9+CB99\n5N5+5x3wWIxddD9pGhJC9D6PwWQMHtx7cfgpSQRCiN6lqp1KBKqq4nS6qLE5sFTbsNocqKraxUH2\nbdI0JIToXXl5UF3t3o6MhIiIJotV19g5XWrhdJmFovIqSiutlFZWU2q2YnM4cKkuVFVFURTCgoNJ\nTYzgxxeMQKNRevBizk6SCIQQvauZuwG7w8XBk0UcyynhWG4pheVmLHYLZpuZKnsVJ4uKMYZqqHZU\no+JCowGn031soDaQEssE0oYlMjA+vIcv6OwjiUAI0bsGDYK//c2dEBITMVfb+GZ/AYdzywkw5VBS\nXUJJdQlWZxWBAVoOfHwlYy7LoDr4KGOHJeF0gtkMJSVQU60jKjiKfqZ+DE+OJyHS0NtXd1aQRCCE\n6F1JSbBgAQB2h5O/rN7C1iOHKKguINgajNEIEXEQEgLfr53A8e3jOb59PPQbQckVn5AYHsmguHiG\nGiOIig9ncFIUEwbHM7RfFIoizUJtIYlACOEzauxO7E4nLtWFQW8gLCgMjUNDTYkGW7GezLUz6sqm\npmi5Y+bVxEWEEh9pYFBiJP3jTOi00gemvSQRCCF8hiE4gJsuGU+UzkqppYYBAwai1WjQaTX878sI\nLKXupp7ISBd/fCqWq6cP6uWI+wZJBEIInzIwPpxzh0YDkJY2DgCXCxbfWV/mvvs0kgS6kNxDCSE6\nxmyu7/bZzf71LzhwwL1tNMKdd7ZcXrSPJAIhRPs995z7EzkkBD7/vOP17N4N8+bBAw/AP//ZZBFV\nhaVL6/d/8xsIlx6hXarFRLBy5UrGjRuHyWTCZDIxbdo01q5d61Vm8eLFJCUlERISwoUXXsj+/fu7\nNWAhhA9YtMj939BQ92RxHbVnD3z8MSxfDqtXN1nkyy+hdsGywEC4996On040rcVE0K9fP5577jl2\n7drFzp07mTVrFldffTV79+4FYNmyZaxYsYIXX3yRjIwMYmNjmT17NmazuUeCF0L0AocD7Hb3tsUC\nBw92vK42TC3xzDP127feCvHxHT+daFqLiWDu3LlceumlpKamMnjwYJYsWYLRaGTbtm2oqsoLL7zA\nwoULmTdvHqNGjeLNN9+ksrKSd955p6fiF0L0tIce8t5fsaLjdbWSCLZvhw0b3NtarbsFSXS9Nj8j\ncDqdvPfee1gsFqZNm0ZWVhYFBQVccskldWWCgoKYPn06W7Zs6ZZghRA+4O2367e1WjB0YvRuK4nA\n827g+uth4MCOn0o0r9VEsHfvXgwGA0FBQdxxxx189NFHjBo1ivz8fADi4uK8ysfGxta9J4Q4S7zw\nQtvLFhbWbw8Y0LnztpAIjh4N4l//qt9veCMiuk6r4wiGDx/Onj17KC8vZ82aNdx0001s3LixxWNa\nG9a9o/bJjw/xxZhq+WpsElf7+FJcE6dOxRYVxaTCQnC5UIGDej2Wc89t+UCzmUmqSu2/8KzrrqO4\no9elqhiee47A7GwCc3LIPXUKcnLq3n7rrYS67ZkzS6muPoav/C/0pd8lwJBOLuTT6h2BXq8nNTWV\nCRMmsHTpUsaPH8/zzz9PQoL7l1RQUOBVvqCggHh5miOEzwr/5BMUh4OgggIUlwsFUICBnn00m2Mw\nsPOrryi88kqcgYEkv/QSk849l3EXXdT+QBQF8/jxFP/oR+T++tfg8QUyJyeA9esj6/bnz5dWhu7U\n7pHFTqcTm81GSkoK8fHxrF+/nkmTJgFgtVrZvHkzy5cvb7GOtLS0jkXbDWozuy/FVMtXY5O42sfn\n4rr00iZfDs7La3uMM2fCPfeg/fOfAdBXVnbJ9dXYHJwus/DwY5U4ne7EMHxcBSUhZXx11IXD6QIg\nNCiApGgjsyamEGEM7vR528rnfpdnlJeXd+r4FhPBww8/zJVXXklycnJdb6BNmzbVjSW49957Wbp0\nKcOHD2fIkCF1vYquv/76TgUlhOhGJSV1myrUNfOgqlBV5R4k1ha33QZnEgGq6n52EBPTrlAs1TaO\n5ZZyNKeEk6fLKa6wcPRwAF/9d1pdmfBJn/P5oWM4XU5U3CuPBWoD6W/qT1CAjiumDm3XOUVjLSaC\ngoICbrzxRvLz8zGZTIwbN45169Yxe/ZsAB588EGqq6u58847KS0tZcqUKaxfv57Q0NAeCV4I0U4N\nxvhUTJhA2K5d9cng5z/H6wltS0aP9t5/9VVYuLBu1+F0UWZ2ryJWZrZitTlwOF04XSp5xZVk5ZVR\nUVVNtauSUmsp+7YmcOTL87CWRVCbnmIHnWL8JZnYbFBe7v7RqkHEhsYSa4ghNbHp1cxE+7SYCF5/\n/fVWK0hPTyc9Pb3LAhJCdKPLLvPaPbJiBZNmzXJ/owf3spHtERQEVisAVf/8mO2X/YyCUgunSy3s\n26eQa84lNdWF1WHF4XIvJ3mqoIq84kosrlJsDjvO/XMp3vMjnDVBDSp3wfBPyNjppH9UDBFBEaTE\nRBJjDGfkwBguGNMfk6HhMaIjZPZRIfzJ7be7p3WorHQ/nDUYyLvpJhKHD4eHH2752J/+FC65xD28\nF3ezjiYmjuDsEwA4Dx7gzU2bsNgsfLz4FhzVtU1MKigqGq0Tjb4GVW9GDT3NrwvXs8bxCwqJa3Qq\nRWclaMR6JowK5cLxlxIfYWBgfDiDkyJJiDLKOsRdTBKBEP7kF79w/9TasYPcu+4isbWHnw4HvP++\n++dXv+LYw4t5d+RMUqZexIjgLWwcMpgt548hJOw43/7lZx5JAEABVcHl0OBy6KHaABXx/JWxDU6i\nEhBWSsw5Gwkb/C3BhDF74kzuu26qrDTWzSQRCCFat2xZ3aYKrI4fw/9ObOOLEf1gxE8JDQWNBk4f\nNnJy97B2VR1BCeqASoZf/l8GpjqJDI7EVj6a4UkR3D43TZJAD5BEIIRo3csv128rCsnDBzEpJwyL\nzYJTdaJRNGgUDSsXzaG+H5JK5MxVmMJd2ApSqTqdgLUsArsllBhbCeNce7iO9xk94jAf3vcYkaEj\nGJQYweCkSCoLTxCk12IIDuiNq/U7kgiEEK3Lza3bVBISuOnScVRZ7Zwus+B0utBqNaxYFojNWv+R\nEhZp5ZafuweFaRQzGuUYeo0eY3AQM9e9yzlvrQSgeNrtPHjdDCKMQXXf/neUnerBixOSCITwB2az\n+6elUf/l5XDddbBxo/uh8Keful8vK3OvFVlrwQIAQoL0DIx3rxBjtcJfGkxXdPwHHTbnTFRAp3Wv\nOxwUoCMoQAdfvFZXLmriWAjruUFhojFJBEL4g+nTYdcu93ZcHDQ1MWR0tPuhMMC6dfWvHzjg1U2U\nRx9tdOjIkd77Dz4IESY9oG86nuXL3T2Yjh6Fc85p37WILieJQAh/8P339dsN5gerM3w47Nvn3q5N\nCABTp9avTZyRATrvj41/P7uHwqwUwAi4Z6X2eLbctIgIdwKQJOATZM1iIfyBZ9NOc81D777rve8x\nSrjO5MmNXrp60WjMZ5JAKsfY89dNHY1S9BJJBEL0dU8/7b3/yitNl2s4ZcRf/9pq1RdfDA5n/cfI\ny/yaFI/2f3F2kEQgRF/XsJ3myiubLxseXr9dUdFitYcOuReWr3UDbzObL6j5+hv2Hy/keH4ZNruz\nAwGLniaJQIi+rrZ9HyCglX7599/v/q/BAL/9rddbNruT4/llfHsgh0+3HGLc+PoP+UCsrOA+905O\nDiv/vZE/fbyJp/+xiT3HPJ5J2Gzuh81VVZ25ItHFJBEI0dfZ7bBpE/TvX/9B35xHH3VPQFdZCX/6\nE6pOhz3MxNGLruQPr/2HP328iVe/2MSCWwOpsdZ/fEwK2koshe5JolUXmYWZ7M7fTV5lPsfzy+rr\nf+YZWLrUnWgaTIAneo/0GhLCH0yfDidOtLl4udnKkSeeJc3pRF9ZwaAN/2bv2GGYk4xsfv52ynPj\n8RxBPPKZdbyxaTYfjpiIaohmQHAU8YZ4xgxM4vwx/esr/tvfzhyiQnFxl12e6BxJBEKIRtZuP8Ks\n11+q21cB1xAj/33ofmzmUDyTQNTkLzhRnIT+4kmMCQonMdI9S+jolFiSY8LqK923z3ua6yee6IlL\nEW0giUAI0Ui4IYjokqK6/fJgIx/+9jFcTq1HKZWk8Qd49uk4YsNTiYsIZUB8OOHNrRHw0EP122Fh\ncPnl3RO8aDdJBEKIRi4ZnYRGdY89sKNjTPV+XHgngb++5GLElBBmjh/YeoUuF/z3v/X7113XpfGK\nzpGHxUL0VWPGuKeT+Pjjdh+qXeSeRsJMKHP5hBySPd5VWf7ySYadc5LY8FAs1bbWK3Q44Ne/do8o\nhsZjG0SvkjsCIfqq2uki5s1z99KprGyxeElFNVl5peSXmHEljGYi8VzFp+zEc9EaFz9+/G12WWzs\n/VqHMcDIiMT+/HpuWv2qYRs2wOuvw4svgsnkfi0gwL3Q/Z//DFlZEBvb9dcrOkwSgRD+IKjpdvvC\ncitH8irYcmI7eaVllFpLqayp5PhRA9s45dUcpGgdXPOnJVRVQXEFKCj0M/Wj2manbu0YrbZ+Oovx\n4+H3v2980pSULr440VmSCITwB+ee67VbZbXz2dbDrN1ylKKaItRglRrVQkQEhMXC1ocfQa1LAiq6\nYAuXPLGc0mId4UHhJEVFEhMSw+iB8VySNqh+FTG9Hmpq3NuffNJ0IhA+RxKBEH3R7t3e+0uWeO0e\nOFHI1kNHySzLxIWLML0Rg8E94Dc/OxDVUT99tEbr5NbnPiNUP4nwuDBSEtyriA3vH924h1BCAhw/\nfuYkB7rhwkR3kEQgRF+0aJH3/vjxXrupiREkR8YwKmIUNc4a4uLi6pebfGI2nuME/vC3XC6cPoPY\n8FBiI0LR67Q0a9Kk+kRQUgKZmRAV1fKCOKLXSSIQwg9FGIO57ydT2bhZwWy1M3r0GHRaDapLyx9v\nr5+PSKtVuO/W/i3U1MANN8CHH7q3nU53N9H9+2HYMFixQsYO+ChJBEL0RZ991moRRVEwBusxButJ\nSXB362y40thf/tLO886dW1u5+wH1/v3u/UOH4PBhSQQ+SsYRCCHqNGzWv+OOdlag1bq7qbpccPPN\n9a/rdI1mMxW+QxKBEAKAq67y3r/ttg5WZDC4/+u54tmsWaCRjxtf1eJv5plnnmHy5MmYTCZiY2OZ\nO3cumZmZXmXmz5+PRqPx+pk2bVq3Bi2E6HqffOK9/3//14nKPv8cysvr9595phOVie7WYiLYtGkT\nd911F1u3bmXDhg3odDouvvhiSktL68ooisLs2bPJz8+v+1m7dm23By6E6DqPP+69P2VKJyt0OiH5\nzLQU8fEwcWInKxTdqcWHxevWrfPaX7VqFSaTiS1btnDFFVcAoKoqAQEBxMqQcSF8w5NPuscNxMS4\n23vasPZwg2EGbN3ayRguvxyys90PHXJzO1mZ6G7tarSrqKjA5XIRUTtxFO47gs2bNxMXF8ewYcNY\nsGABhYWFXR6oEKKNVq50r0qWmwsvvdRq8Y0bjahq/X5SUhfGMmIEXHRRF1YoukO7uo/ec889TJgw\ngalTp9a9dtlll3HNNdeQkpJCVlYWixYtYtasWezcuZOA1tZHFUJ0vZKSdhVfuHCI1/7Bg10ZjDgb\nKKrq+V2geffddx/vv/8+mzdvZuDAgc2Wy8vLY8CAAaxevZp58+bVvV7u8eDoyJEjHY9YCNGiSZMn\n140LdikK3337bbNlT54M4JprxlA7kjg42MHXX+9utrzwTUOG1CdzU+2Mr+3QpjuC3/3ud7z//vt8\n9dVXLSYBgISEBJKTkzl69Gi7gxFCdC1naGiL799880g8p5N49dX93R6T8D2tJoJ77rmHNWvW8NVX\nXzF06NBWKywsLCQnJ4eEhIRmy6SlpTX7Xk/bsWMH4Fsx1fLV2CSu9unNuPTjxjV5XlVVKSm3YTbX\nzxuk0aq4YmBfkQ6H04VGUYgwBjEgLpwJQ+LrZxjtAfK7bB/PFpeOaDER3Hnnnbz99tt8/PHHmEwm\n8vPzATAajYSGhmKxWEhPT+faa68lPj6e48ePs3DhQuLi4ryahYQQPUijqV8T4Mw00A6nixP5ZRzL\nLSW3uJKCUjPP/HYCEHjmIJWJ1/yXf353CJfqwqW6UFAI0YfQ39SfcMMUUhMjmjydOPu1mAheeukl\nFEXhogZP/RcvXszjjz+OVqtl3759rFq1irKyMhISEpg1axYffPABoa3ckgohuonTWbeZlVfKtxv2\ncSSnmGJLGSXVJVTUVFBitlBedJ7HQSoTr95KTQ2UlUFFGeg07rUHQgKCMIZIx4++rMVE4Kr9VtGM\noKCgRmMNhBC+4f2vMtl59CQnyk9QVFVEUIiD8HCIS4BPb34Ez2cDYeM2kPGdlYTwKCKDI+kfF0F0\naCSjB8Yxa2JK43UHRJ8is48K0QdZbQ4yj5/mYNFBKmoqUFExm8FiVsjNmITTpvcq/6s7qgnRzybG\nZCA10b3wTEp8OIEB8hHhD+S3LEQfFBSgY3baIOKOGThdZqlr89dpNTz+xHQ87wauWfANv/vxDGJM\nIfLB76fkty5EH3X+mP6cP6Y/qqridKloFIWRI717/oSFOXj4thCSY8J6KUrhC2ReWCH6Eq3WPfd/\nbCxkZADuaWB0Wg1ZWQqHDnkX//LL73shSOFrJBEI0VdUVLi7jTqdUFjoXi/YQ8NhQO1edEb0WZII\nhOgrHn3Ue3/+/LrNBQvqhxaA+8ahDZOSCj8hiUCIvqKZdYprauCVV7xfa9hEJPybJAIhznIul0qZ\n2YrLY95/Ffjf3pN8eyCHuDinV/kRI2DQoB4OUvg06TUkxFnEanNwIr+M/BIzp8ssnC618OXWUoYN\nhaU2W105l6Lw9w0bObknmfLyRI8aVN779DTfHXa3Ex0/bSY6LBDh3yQRCOHjys1W9mad5nB2MSdP\nl1FSXUplTSVfrxnDD1sm4x4ToPIyxQznIKPZR4kawSePz8BREYPnmIHpP9/OS/85gUt1J4Ky4jJC\ndaEERfZjysjkXrpC0dskEQjho2x2Jxt2ZbEl8wT55gKKqooorynHYHSy9a+/oPBoKvUf8gplRLKN\naWxjmvulCu/6NHobMdPWUQhUVILNBuZKC/1Ck8kvMffglQlfI4lACB91KLuI/+46wJ6CPdhddhQF\nTCb4fPFvsZZHUJ8E2kLlmheewWaDygoNxoAw+kVE4dQ6GZYYyczxA7vpKsTZQBKBED4qMcpIREgY\n48KH8ZvlD9H/1HGCdNU47Xo8m3vQ2AlNycSaOwinNRTUxn1AjEO/ozI/hsHx8YT3C6d/TDiDkyKx\nl+cSaQiUSeX8nCQCIXxNfDwUFBD11ls8PGo02kkzUYDXmN8oCYQYrcx/5jPsLjtwAI2iQaNoCNQG\nEqQLIjggiG37T3HLZeOJDb+C2IhQBsSFExLknnRux47iXrpI4UskEQjhSx56CAoK3Ns33YR28mRU\n4Gke5TGWeBRUGTSsiidf/AGbYwgRxiCSosPQaTXotBqMwQFEGIMxhgSw6fsT0vQjWiSJQAhfUVKC\n+txzno0+bLMY+T2b2Yr3IjIDJhwg7fqv+XS3exWxl+6/Fb3DjqLRwMUXw+ef15WWJCBaI4lAiF5W\nZray81AuF6YN8hrhuYVzOH//F3g+FJ7FF5wzdh0HbzZQUAaqCtEh0e4kAO55JL75poevQJztJBEI\n0Ys27j7Ohl3HuPeea7z6AK3mGn7GGjyTwM95hzeYj7LXwU+tL5AaHkF0SDQDY6O9Kx0woEdiF32H\nJAIheond4WTj7uPszP2OZ+dcyrL330cFlvAoj/MUnklgDv9mFTeioLJ/3DR+NPYCBiVGMrRfFHEh\nOu+OpL/4RQ9fiTjbSSIQopfodVrShiWi0ZyDJXY0s/UPkvuPYA4w0qOUikbnYuaf4fuLcoiLCGVo\neCij9dr6Is89513x73/fI/GLvkMSgej77rgDXn7ZvR0SAhZL78YDsHIl3Hknl50zmElDEpg2Rc/x\nY4F4DxJTGTyikqdXZpMUPYah/aIwBAc0ruuNN7z3A2XuINE+Mvuo6Ptee61+u6oKVqzonTjeew+X\nRoOqKKh33cXKj77l6be/JjE2tMkkcM4vPmDUTa/xwY5NvPvNt6zPONZ0vb/5jXvIsUYDen3TZYRo\ngdwRiL7PY1ZOwN10ct99PXb6/nfdTcz2bah4f/O67ObLGVxZSMMEACrX/fVJysqgzArhQeEEagOb\nvhsAuOsu948QHSR3BMI/3Xtvj5ym+PKridm+DQXvj3s7Wn5S+TkNk4A+7DSzFj+NoyKK/qFDODfp\nXC4eNo2bLz6H2Wmp9UWtVvjggx65BtH3yR2B8E8vvwwvvNDtp6nxWCwG3N/3C4lgAKewEuL1zpR5\nO5k+pxBD4HQGxIUzKDGCwUmRxEaEoihnEsbPfw7vvedxmNrt1yD6PkkEwj/V1PTIaezr1vPYhxtI\n/+1P0LpcnD/tO7ZuGU/DO4HfPv0dt18/gNiIkUQag9FomplZ9MAB7/1vv4Vzzumu8IWfkEQgRDca\nEB/O1ef2Z9fWbTzwwAS2fq2l4fTRZrNCaOiktlW4bRsEB9fvz54N5eVdF7DwSy0+I3jmmWeYPHky\nJpOJ2NhY5s6dS2ZmZqNyixcvJikpiZCQEC688EL279/fbQEL0W6q2mtNKFabg/zSai6/ciRfN0oC\nLv784Xbe+iKDDzft58CJwtYrDGowXXRFBSiK+6fhe0K0UYuJYNOmTdx1111s3bqVDRs2oNPpuPji\niyktLa0rs2zZMlasWMGLL75IRkYGsbGxzJ49G7NZVjwSPuanP4WwMHd30m5IDKqqUlRexfb9p3j3\ny708v2YrT7+9iasvP4fiwhC81hDAxZwly1m7fyP/ztzEmoyvWbMpE6vN0fqJmlt5voeau0Tf02LT\n0Lp167z2V61ahclkYsuWLVxxxRWoqsoLL7zAwoULmTdvHgBvvvkmsbGxvPPOOyxYsKD7IheivTwf\nsnaUVuv+9u2o/8CusTnYtv8U3x8rILe0lKW3zSE6LoFrZj3H9+/ecKZUfRLQ6Oxc//JSysqgshKM\nAUaiQ6KJMoYQ6DliuDk7dkBEROevRYgz2vWMoKKiApfLRcSZP8KsrCwKCgq45JJL6soEBQUxffp0\ntmzZIolA9C0aTf2dREAA2GwUl1fx2n92cawom1MVp1jxgnsa6cKCfhS8exENHwprgyu48OEXqSqK\nIj44kpFJkSRGRDB5eBLnjEiq7x3UkvBwiI2F88+Hf/6z/vWAZsYZCNGKdiWCe+65hwkTJjB16lQA\n8vPzAYiLi/MqFxsbS26DbnNCnNU++cS7Ocluh6AgcjJPkF9eQk5FDlaHlYRCCw/zLMu5Hxee3+5V\nYvoX8asdF1v3AAAgAElEQVTHdxEWPIPEKCODkyIZlBjh3T20rWoXr/E8Lj6+w5cn/FubE8F9993H\nli1b2Lx5c5v+aFsqs2PHjraetsf4Yky1fDW2syGuiZMn122XnHsux198sWOVJiYyLCUFQ1ZWfSNP\nTQ0jRg0gZuUaLJpEItOPcA47OMgIjwPdI4Vn/eZ1Zk+K4ZzB0RiD9SiKHawFZP9QQHbHIgJgEvX3\nHPlpaZzq4O/kbPhd+hJfi2vIkCGdOr5NI4t/97vfsXr1ajZs2MDAgQPrXo8/8w2koPbbyRkFBQV1\n7wnRmxTcf+QaIHznzrrXkx57jFGXXdbq8YNvuYXBt9wCwKH336cqMZHa+wIF0NXUcPivR3h14bU8\ny6JGSSA4KofpDz6B1lTAvlOnySutbv+3/xaUTZ2KLTwcl0bDqV/9qsvqFf5FUdWWu0/cc889rFmz\nhq+++ophw4Z5vaeqKklJSdx9990sXLgQAKvVSlxcHMuXL+e2226rK1vu0dfZZDJ15TV0Sm1mT0tL\n6+VIGvPV2M6auMrL3e3pta68EjZtcj+hrdVa76EzH9rqmZ+lb23kjrt/TGR5CQAFmEigFM9nAQYq\niaKIYUvfxG6HynILUUFRpA2axC8vSyM10Xce9J41v0sf4atxdfbztcU7gjvvvJM33niDf/zjH5hM\nJvLz88nPz8dyZhpfRVG49957WbZsGR999BH79u1j/vz5GI1Grr/++nYHI0SXuvxy7/2334YZM7xf\nO3WqTVUpuBPBxqxvuHnR3ZSHBGEN0DVKAhPZyVrmMOKp94kKSGRk9GjGR43jgoEjuPvqqT6VBISo\n1eIzgpdeeglFUbjooou8Xl+8eDGPP/44AA8++CDV1dXceeedlJaWMmXKFNavX09oaGj3RS1EW2zf\n7r1vMsGnn3o/YJ05E44ebfr4HTtQqf+Y352UhMPlwGyGeQ88zMYnHqNhr6CJfzjF54EPMDsultTE\nCAYlRlCef5ygAC1xkYYuuzQhulKLicDlcrWpkvT0dNLT07skICG6jNPZepljzczxD3DppV5DwD5d\n8jbTVRcKCk8vuICGSeCxv+7g1nkziI0IJTiwfl2AHSWdeRwsRPeTuYZE3xUZCSUlHT/e41gFSL95\nBk6XytQpCk6H9wPfJUsUHr1jMkKcjSQRiL6ruLhLq1MUhZdfUmjYc3D4cHj00S49lRA9ShKB8D9P\nPw2pqfCzn7VcrrZH0YsvwrFjVFbC3Xd7F9FqG88MLcTZRhKB8D+PPNK+8meWgQxTXHh3tHOxeU8u\nWzNVggP1JEQaiAkPbX4tASF8lCQC4fdqbA4Ky6soKq+itLKaMrOVvBIz+cVmysxWyi1W3npkLt5J\nQOW2P7/Hq1+4O1ToNXrCAsM4b2Qqc88b1uR5hPBVkgiEX1FVleKKao7llPBDXin5JWZKKqtYcd85\nVJXH0nDRGG/1fYjOv+cVyrS5VFZCdTUoKAwIH0BRuYyoF2cfSQTCt5x3HmzZ4t7uzJoBOp1X91Fr\njZ2tmdnsOVZAXlkZJdUllFaXUmmrZPvbs6kqN9ByEqilEjVqB1EDcyku0mAMCCMhIpLY0FhGD0jk\n0snNrBUghA+TRCB8RxfOweOZBFTghQ+2klV8ipzKHKqdZkwmeOidfzBx7zGmuC7iuzYmAV1wFVfd\ndpDQgLGE9zPRLyacQYmRjBwYQ7wMGBNnKUkEwndccQX8+9/1+/37w8mTna5WBfIrit1JwF6NioY1\n96SjBWoI4DipXqXDKeZVFhBDEVEU8WDkMpJ/bSI0sZrzx/QnNmIGseGhJMeEERKkb+asQpw9JBEI\n3/HZZ953BdldMyJXBS4cPZwhhYmUW2pwqSpVoWEYLRW8zY2UEFVX8sEnD/P0k6PRnVmBTAXW5FxG\ncKCuS2cNFcKXSCIQfZ42KIirzh/u/aJpFc6rruYPPODxosKyx4bB4w6PV5Bv/aLPa9N6BEJ0O0Vp\n+hnB6NHtr6thc1JTM+HOncsafsIh3AlCi4Nf/rL9pxKiL5BEIHzb00+3/xiPudkB+Pvfm66aRXXb\n57Kt6WIyi67wA9I0JHxTZ7qOjhnT6vE6HTgZ497GzkM8B5zvvWgNwJtvdjwOIc4SkgiEX/KcodqB\nnrnqJ+4do7FzSUiIs5A0DQnfExjYrdVHR3vvG6T7v/BzkghE72u4ZsBVV3Xr6RrOTt2wNUgIfyOJ\nQPS+khLQePwprl7t/f6UKbBoEV3hvPO892VogBCSCIQvGDzY3Wivqt7t87VdSrdvb1/vodrjFAWe\nfdbrrdppjGqVlnYibiH6CEkEok+Je+YZ7xeuuKJuc/HixuVN993q7kLkmTwUpfGDBCH6MEkEwnc1\n7MN///2tHhL/2WfeL4wZU7f5xBOeb6hs+9aO67XXml7kvouXuRTCl0kiEL7LbPbe/+MfWz1EY7PV\nbatAhaWGYzklvPxW9plX6t/9ZM9GnBr5JyCE/CsQfUrDZ7/PvreJy64u4Y6bk/BcWKb/hf9he97/\n2NU/uemKBg7sviCF8DGSCETv82ybj431fi8goE1VlFZW8+WevLqPeht6nuAxnvnVbI7uGETDP/Ur\n78ggKQmW3HwNDYePqQB79nTgQoQ4O8nIYuFbdA3+JPPzISoKMjJg0iSvt6pr7GSfriD7dDn/9++d\nbMnM5FO+4DRJnGQANQQ3cQIV4r7n++8CSY2NIyUuAnjeq4QC7hHGQvgJSQSid/3hD977b7zhvR8R\nAS73AvGqqpJfYubAiSKO5ZaQXVhGmbWcUksl7z95HTiv5VSzK42p7p/wY1zyi7389ILLSIo2kpIQ\nAb/p6osS4uwiiUD0rqee8t6/5JImi+UWVbJ+xzEOZOdRWFVISXUJlbYKDEYXnzywCNDS+AnBmUYf\nrZXLF71BQlg0OoK555rLGDEgxqOYzC0k/JskAtG72jC/g6qqrP5qHztO7iO3MhcVleBgiI6B93/z\nOO4EUP8g2M3JzxZ+yIBUE3EmE4MSZzG8fzSDEiPQauXRmBCeWv0X8fXXXzN37lySk5PRaDS82WBa\n3vnz56PRaLx+pk2b1m0BC/8UFRZCojGRoVFD+dHR04RYB59JAhq8k4CDt9fvZX3GCR7+5Qzu/8l0\n7r12CnPPG8bQflGSBIRoQqt3BBaLhbFjx3LzzTdz0003NVq3VVEUZs+ezapVq+peC2hjTw8hePBB\n+NOfoKam2Yl/FEXh+ovHYE+ZS2D2SZxo+Ae34v09RiUiwsozL33NDbMv7ZHQhegrWk0Ec+bMYc6c\nOYD7239DqqoSEBBAbMNuf0K0xbJl7p9W6E4cR5t9kmqC+TnvsgvPHkQqYycXsysjBuVnoP4MXIBW\n2v6FaJNOPyNQFIXNmzcTFxdHeHg4M2bM4OmnnyYmJqb1g4VohtPp4lRhBVn5ZRSUmDldZiFbk856\n19V8z/i6ckM5iP7ybIbMOISSUd9IJJOKCtF2nU4El112Gddccw0pKSlkZWWxaNEiZs2axc6dO6WJ\nSLSL0+ki83ghmcdP80NeKSWWckqtpVTWVPJDZiDfudLx/Ii/nz/wLA/BWpVLpqR71eXSaGS0pBBt\npKhq2++fjUYjK1eu5Kabbmq2TF5eHgMGDGD16tXMmzev7vVyjwXFjxw50sFwRV+VW1LFpswCciqK\nKbYWU2GvQBdoIyzMQUiIk3/+9kXqnwmoLGUhC1mGCvwQHsn9dz/OP5+6ty5NlI4YwbG33uqdixGi\nhw0ZMqRu22Qytfv4Lv/SlJCQQHJyMkePHu3qqkUf9v3xUg4WH+dIxRFKbCWgtQNQVqZj9+4gGjb2\nXNr/Q1QgNzyWd597jWX/We3Vd+jYX/7Sk+ELcVbr8nEEhYWF5OTkkJCQ0GyZtLS0rj5th+3YsQPw\nrZhq+WpsXRaXTlc3BfQk4Nm3v6aiugqb0z2DqEbRoNVoWfbsxXhNGJdSRdCGLRRHhJJgDGaRRoG7\nrqurVgHSLrqoc7F1IV/9PYLvxiZxtY9ni0tHtKn7aG1Tjsvl4sSJE+zevZuoqCgiIyNJT0/n2muv\nJT4+nuPHj7Nw4ULi4uK8moWEaJLHOgAK8PD151NmtlJd40BRQKfVoNdpWXab1uMghRM/hAIN1io4\nfZqSGTMI27sXnfQWEqJdWm0aysjIYOLEiUycOBGr1Up6ejoTJ04kPT0drVbLvn37uOqqqxg2bBjz\n589nxIgRbN26ldCGi4oI0QpFUYgwBpMYbSQhykhMeCgp/bybhfT6Zg42mfjh1VfZvX173dxEQoi2\nafWOYObMmbha+Ie1bt26Lg1ICE9lZd77HuvOCCG6iPSwEz7LY7lhIUQ3kkQgfNbatd77soywEN1D\nEoFonzFjmDR5MpMmT252bqA2e/11iIx0b193nddbH37YuHhtUSFE15JEINpn3z6vSZ+bVFwMBkPr\nX+Hnz3eXUVVYvdrrrWuv9S763HMdiFUI0SayHoFouyYmHWzE8y4hORmqq9t9mpKSxq898EC7qxFC\ntJHcEYi281iLok099a3WDp0mKsp7f+zYDlUjhGgjSQSi43po4Nb33/fIaYTwW9I0JDrEATQ5tis4\n2Ls56D//gTPrWQBUWe0UV1RRsno1paPSKNUHU2N3YrM7sTmc3PWTUUAAtU8hAgO78SKEEIAkAtEe\nZ+4Aqk0mMt97jyZnW8nOhujoul3H5Zez4dujnC6zUFBq5tn7hlGUG0ah+hBjcY8WcwHXPvMCTtVJ\nTfUEPOcV2vhdDt/scWJ3uHC6XOh1WmJMIfSPM2EMkSwhRFeQRCDaLfPLL5t9rzQghHCPfQV4e/Mm\nLDYL/3z4DkALKMRQSn9OMIxDFBDDnoXjaPzkwcWrX2zEpbpwupyoqGgUDYYAA4Nikrhr3jnodVqE\nEJ0jiUB0CavNwZc7f2D7wWx+G2okL8zEvb++DYPRRUTESTY89ktqk0CtkwzgJAOaqVHlwieeojrY\nPTddRQVUVYGCQqIxEYs1ttH62UKIjpFEILpEZtZpvvj+APtO7+Pr++8DIMzgQlVh//cKZaf64dnk\no8WJs8U/PxX0FsrLQ6m2aIkIjiA5KpLokGhG9o/n4kmp6LTS10GIriCJQHSJ5JgwIkPCGR8/HqvD\nWteMo1W0fPrCj/C8E+jPEY4wmh9IJZORXMsaGnZgO+eGjxhsmIAhwEBErImB8eGkJkQwYkAMkWHB\nPXtxQvRxkghEy/71L7j66vr9zMwmi8VFGnjkhgsoLK+i3GxFBfRaDYOTTXiPQ1Y5oBuF3uFgOIcY\nxiHuWvkxQboggnRBhIUEUVhexYyxY4iNCCU2PJTEaKM8CxCiG0kiEC3zTAKt0Go1xEcaiI80ADCg\nUfO/yp//r5A952XB0aOMuudWAvPzuOvKGYQbgogwBhOglw98IXqaJALRPiNHwpnl+lqyeTOcOunE\n/YAYNDiJidNy922xZ+pJhrnHABjWXbEKIdpEnraJ5u3f3+FDL7gAXGeSgIKLR1lCfn5XBSaE6Epy\nRyCaN2qU9/7Gja0eoqoqGo0L93cM97OBx3mSdJ7grc+vxu5wotEoRBqDGRAfzvjB8V0ethCifSQR\niLabMaPRS6WV1RzPLyO/xExBqYU7fzoUCKY2CVzEFzzGUwB89MWHuCIjUBSFYF0wKREpdQ+DhRC9\nRxKBaJ6quqcC9ZgX2uVSOVVcxYnTZrae3E5eaTml1aWYbWaycsxUVYzDc7zAKm5Ei3vN63dWPM2c\nxx7jbytX8tov78KQMIKwUJkmQojeJolAtOzM4jKqqpKZdZoN32WxPfMQpbZSNCEarC4L4eEqxij4\n3+8fxzMJgJN4CuqqCnS5+Em2haFFRTz73GJgMcrnn8Mll/TwRQkhPEkiEG2SlVfG219+x8Gig2SX\nZwOQHGUkJsQ9BcSWD0fRcN2yX77wLuq97m0VqDEYufXjN7xLSRIQotdJIhBtEqDXotfqSDAk4DA7\nAIjSxKKz6dAqWg5/eg2edwNDR1dwy6wZ/HCqmPhIA4bgAIKh8+scCyG6nCQC0SbJMWHcddUU8krM\nbNuxG4DRo0YSqNdy5y0heN8NKBzaawJMvRGqEKKdJBGIxmq/tW/c6NVTKC7SQFykAUeZe72BtDH9\nAcjY5n340qU9EaQQoqvIgDLRvJkzQdPyn0jD9YUBFi7snnCEEN1D7giEN32DBSgXLGixuEfPUqDZ\nOenqFRXB4cPw3HMQEND++IQQXa7VO4Kvv/6auXPnkpycjEaj4c0332xUZvHixSQlJRESEsKFF17I\n/k5MTSB6mcPhvf/yy80W1TYxP9zIkc0UVhT3T3Q0TJsGH30Eq1d3PE4hRJdpNRFYLBbGjh3Ln/70\nJ4KDgxutCrVs2TJWrFjBiy++SEZGBrGxscyePRuz2dxtQQvf4HJ576sNV5oUQpwVWk0Ec+bMYcmS\nJVxzzTVoGrQXq6rKCy+8wMKFC5k3bx6jRo3izTffpLKyknfeeafbghbd5KmnvPeDm18AZvLkCV77\nrTxKEEL4sE79883KyqKgoIBLPAYFBQUFMX36dLZs2dLp4EQP27XLe7+qqsli7pY/7z8dp7OVuh95\nxHs/Lq5doQkhuk+nHhbnn5lXOK7BP+rY2Fhyc3M7U7XoDf/8Z/326dPNFrv55ol4Dh4LMVh5/T8H\nsdocWG0OnC4VvVZDTHgok4YmMKx/NDz9tHe/0hbqF0L0rG7rNdTwWYKnHW1Y2KSn+WJMtXottpMn\nqbE7ySutpqjCSlFlDZ99EANM9So25a7lvLHRiVN14FSduFQVraLFqDewfXcqN85IRVEUJuE9E9HO\nbrouX/1d+mpc4LuxSVxtM2TIkE4d36lEEB/vnku+oKCA5OTkutcLCgrq3hNnn4oqO1kFlZwsspBb\naqHSZsbisFDlrGL7+kfx/Dg3pXxPRFIeqqpgNmuxWLTYq3WE6Y3EBccxMNZQ96WgaPJkojMyAHDI\nQwUhfEanEkFKSgrx8fGsX7+eSZMmAWC1Wtm8eTPLly9v9ri0tLTOnLZL1WZ2X4qpVk/Hll9i5r87\njnEgu4SiqiKKqoood5RjMDkJD4PPfvUYDReiH3/rPzCbY3Ha9EQGRzIoIYqokEgGJ0QzfdwAUhIi\n6ot/+637vyUl6CMj6eqr8tXfpa/GBb4bm8TVPuXl5Z06vtVEYLFYOHLkCAAul4sTJ06we/duoqKi\n6NevH/feey9Lly5l+PDhDBkyhCVLlmA0Grn++us7FZjoWXnFlfz1k+0cK/mBvMo8VFQiIiApxv3+\nqlsfxXPVMVAZ96P/MsKUhiHAQERoGIMTIxmSHMngpEiMIS2sMxAZ2c1XI4Roj1YTQUZGBrNmzQLc\n7f7p6emkp6czf/58XnvtNR588EGqq6u58847KS0tZcqUKaxfv57Q0NBuD160n8ulYrHaqLDUYLU5\nUIEIQxDx0WE8ibvt/sNzz2fNNT9Ha9eilun44Mmrcf+p1CeBwNAynnk0lYQoI/GRBmJMIWi10twj\nxNmo1UQwc+ZMXA1HDjVQmxyEb6msqiGnqJJThRXkFFZQUllNRVUNVnsNNqcNh8s9injhIwu8GnzG\n6wNh8gwCdFpW/S0Ya3kk3o95HWzeeMznbo+FEB0jcw2dTWbMgK+/rt9vMJS3zGwlK6+UrLwyjueX\nUVxppqKmgoqaCiptlVTZq6hx1KDTqwQGuqcVKi2F2LIir3oK/voaETYH+fkqH77dn4ZTTGdkfN99\n1yiE6HGSCHxdXl6zi7moqkr26Qr2/FDAsZwSCsorKK0uZeNnERz8/OIzpVwYEnK4+c+vExoKOh2U\nl8OpU5CdDVWVehrW/tqGjThcDlbdN5+Gq46pKnz7rYpLVVFVtcVuwkKIs4MkAh83ce7cJl9Xgc9/\n+mu+vHguBZYCiqqKsKlVVFfDwc/Tqf8A12DOG8DKn6Q3ONoFIUWEz/mzV71OwBr8A+/e/jjeo4dV\nFq7cztK37WSdOImKyldHa4gKC+GKKUNIignrqksWQvQwSQQ+riomBkNhYd2+euan0BDGyvFDqSnK\nIC4ORgwGgwGenOW5gLwnpcG2BqriqfjwKc7jSoZzkMEcYSkPYbndSMNuoiPmvs93JVnYC+2UVZQB\ncNx5nAGmAUSGBXPtjOamHRVC+DpJBM1QVZXKKhvVNXY0GoXQoABCgvStH9jFDq5dS9rkyai4v60v\neOhZTobYiYhykJJQSUSEexbQmhp4ctYiGn6A12u6CceFlq1MYyvTmolAJXLgD1x26wGsVnezUo3d\nhbMmiIHhAxkY0Y+xqTJvkBBnM0kEDTidLjKPF7I1M5sThSXYnDYUFAJ1gYzsH8fUUckMjA/v2bZx\nVeXt9d+zds828sx5hAa7lw344Qd3Avjlmo+J+97AH/BsElIxpq1m1q+2Unl0ABuW3wVoz7zf8O6g\n2RODYmPW/avYtw90aghRwVH01ycSGRbGlVPSuGBMf0KDZYEZIc5mkgg8HDpZxL+3HeFEcT7ZFdlU\n2EoJDFJRVbDVaDhWEseuH5IZHB/HNTNGEm0K6bHYBidFMujUIBKMCdQ4atBqtOjLilm+eAFONMxh\nHV5dPBUrk648hKUglqBwJ5c/sZL8AgWrs4oqtYQaVxVaJY5Tb9wH9trF5xsmBSc3LF+FUR1BakwY\ncaZwhvWLwl6eR0JkCOeeM7jHrl8I0X0kEZzx/dF8PvhmH3vy92LXltOvP4yOq59n32ZzkZubx/c5\n+ZyqiKfMYuUXs8eRGG3s2kAUBe67D/74R6+Xp4xMZkxqHKWV1VRW2dDrNIQE6mHxAp7icb5gtkdp\nF4+8shmdZjrBAXpCg/WUVVr5znWcIruZCMWAXheGXqtjbPrL6LV6ArQBhOpDCQ0IJUQfQmGJjWmj\n+pEYNZ2kmDCSoo1Em0JQFIUdOyq69pqFEL1KEgGw90Qph3YXsadgDzGJVQwc2LjHZkAADBwI/fqp\nZGbm8e1JG/b/OLnh4rHe8+l0Ru1JV6xw//zrX5CYeOYtBUNwAIYGzTALeJlXua1u30QJz689QYXF\niLnaRrXNTo3DgTEogIEx8cyMH05lVQ2KAoF6HSGBekKC9IQG6YkJDyU2PJSY8BAijcEyUlgIP+H3\nieBIbgVf7c+mgAKSB9bQr1/L5bVaGD0aDh4sZkfOLpzrXdx59bnd00x0Jgk0pKoq+SVmBiTqsbOA\n2iadi/iCz5jDjXufx+qwYnVYcbrcK8YE6YIYGjWU6WPHMmFIgrsijQaeeAIee6zrYxdCnDX8+itf\nZVUN/ztYQFZlFv1SWk8CtTQaGDECAsLKOFR0hI++OYDL1ckFe/Ut90iyO5wczi7m0y2HWLFmKz+e\nfxS7LZDaJJBIDu9wPYE4uOTApwwbY+bcqQ7OnaISEami1+oxBIZiCg1yV/j3v7tHhz3+uPtO5Mc/\n7lz8Qoizlt/eEaiqymdbD3O8Ige9oYKkpPY17ygKDB0KGRl57Ms+yZbMGM4f07/jATkcDQNEzcgg\np7iKYxv2cSSnmCJLKUVVRRRXFbPtswfxfDj8NtcTQyEuIG9cCori7lVUeFpLkiGZYf0H8ZMZo0hN\nPHOdv/qV9/lefrnjsQshzmp+mwj2HCtg57ETFNrzGDHY2twsDi3S6WDYMJVD+w/x350mhiZHERvR\ngVlXExK8dp3AdwdzWLPlBLnlZThDDlJYVUhwqIOoKFj3+8brAqT2zyT99XSKiyEnB6p3BZJoTGRK\nUiITBiUxa2IKEcbmF6MnNrb9cQsh+gS/TAR2h5PPM45xsOggyclWAgI63qwTGQkRMVaOFh9j/Y5o\nbpw9tv2VfPopnBk0BrB8+WpOfrmFfbmZ1CgVjBgQxuR4CAwE91LQ3lM/pD3yAL9PuIvS/ykY9SYS\njYkkpcQzcUgC5wxPIqoHu7kKIc4+fpkIdh3J52RpHkqgmchIe6frS02F7dvyyDyRR15xCglR7exS\nmpZGpcXK+19lMuuWK9lkPY5OqSKm32nCwx0kJ9fP4/PKDZ53AypQw0BTKlFEMSwxin4xEUwYHM/4\nwfEEBjTz673hBu/9pKT2xSuE6FP8LhG4XCpbMrM5UXaC/kPA3vk8gF4P8Qkusiuy+d++7HbPu+Ny\nqfzfZzvZk3OI9357OyNGVBEeDidOODGbteTm4p5M7iA0fL6/+I2tDEq8gKHJUQzrF4XJENT6CT/9\n1Hv/1Kl2xSuE6Fv8LhEcPlXMiaICHNpKoqPdszx3hX79IGN7Pnt+yOPSyYNaXqqxAZvDSUWVleLq\nYpxOhX37VFSXgtUSSaA2ECU6gWB9MDuW34Dn3YBOV8PDP7+g6W/+ycnuhwVqE81eFWcGhP3973DP\nPe2+ViFE3+J3ieDbAznkVOaQmNjsNP8dEhgIEVEO8irz2Xk4j5njB7b52KAAHZdMGkyMyUBJZTV2\nhwO9VkdpUQFhIXomjhlBdWUwr3vdDSjY7U18+7/oItiwwaNYg4v0TAy33ur+EUL4Nb9KBOZqG4dz\nCimuLmJIfNfXn5gIxw7ksfeHgtYTwY4dMHmye7hyVhbnjenPeWP643Kp2B1OAvRadu7cCUBa2vBG\nn+fNDjv44x9hwoTOXooQwo/41YCyI6eKKakuwRTuam38VoeEh4PVZSa3pIyi8qqWC0+e7P7v8eNe\n39o1GoXAAJ3X7KbunkLebLZm6h0/vn1BCyH8nl/dERzOLqa4qpjobuokoygQFQXF1cUczi5uftqJ\ndk7p0LBTT6tJTFXhb3+DV15xT5JkNLrbru6+u13nFUL4B79JBA6ni6O5JZRUl5AS2X3niYqCgqxi\nDp8qZtroZuasWLLEez+k+X7+7bob8HT77e4fIYRohd8kguP5ZRSZSwkMsRPUhh6WHRUZCYcOlpGV\nX0J1jZ3gwAZf35v6Om+xNFvfVVd5t/fr9TYys8ooKLVQXFFFuaUGS7UNjUYhLCSQof2imDw8CZ3M\nHCqEaCO/SQRHThVTXF1MVFT3nkenA2OYk5KqUo7lljI6pcHUDQ3nFHrttWbramoU8YMvf81L/ynD\nbKN4+VYAACAASURBVDNTZa/C6rBic9rQKBqCdEEMOD4ARVGYMjK5qy5JCNHH+U0iyCmqpNxazsDw\n7j9XeDiUl5aTW1TZOBGoqneXzltuaXS8udpGZtbpM3cD3qOIt576H2EmlTATRAZDcLD7MYDdDrt2\nVREZHImlui1tR0II4eYXicDlUikoNWO2mTEYuv98BgPkFpjJLzE3XUBV3YO9PJ4CO50u9p8oZPfR\nfA7nFPH9kULgp16H3fvxsxiN9aum1bJY4MgRiA9NYFhsCpOHy5QRQoi284tEUFpZTUW1BV2As1u6\njTZkMIDZ5k4Eqqo2vdD9mSRQZbWz83Au3x7M4WRxAXnmPEqqi/nymUfwuhvQ1mAyeVfhcMCJE5Cf\nq2OgKYWRCancfOl4wkLbPqpZCCE6nQgWL17Mk08+6fVafHw8uU11d+kl+SU9dzcA7p6aLsVGWZUF\nc7WtyekmSiur2bz3JLuO5pFTnsepilMQYCExEULs0HCIR/oXz9Zt22zu5wc5Ocr/t3fv8VHVd8LH\nP+fM/T6TZJLJjSRAIBAQgUghVUCrVNs+bq1aq92qaHVrV6uo3eqju9XW6qP7PHbbrbdit0tfXVfW\n19pnL3ar7goKj7ItAnLVikRIgAy5zUxmMtdzfs8fgZAJCdchZzS/9+uV1ytz5pyZ73yT1/nO7/c7\n5/ejzBaipbKe1pn1XDx/8rGD05IkSSdQkBZBU1MTa9euHXpsMpkK8bIFE+5LkMgmcPtOvG8hKMpg\nqyCRSdC57QM8C+cMTe2QzuR48729vL1zLx/37WN/bD8ef5YpMyAQGDz24QtHzDCqpPi4M0KZ209H\nB3QfMhF0lnNusIaZtVUsa5lCddA7ZjySJEnHU5BCYDKZKC/ihU2OtAjKT2PNmNPlckF8IM6UhXMG\nNygKAw4Hz/z9f/LH8F7aIm0EyjLMbcm/jeA3j1/EyNbALb9+nPBBPwfiVqo91SysqaK5LsTCmTXU\nh/yjdz1JkiSdpIIUgj179lBdXY3NZuMzn/kMjz76KA0NDYV46YLo60+SzCaPd99WwTmd8D//4q68\ndcQcySTr9vwezRJh9rmDN/yOtPV355PXGnB10NcRpNZdQagsyPxpVXKxGUmSCuqMC8HChQtZtWoV\nTU1NhMNhHnnkEVpbW9mxYwclJWfxFt5TkEhlyWgZrNYxdti3j7tv/jue/I8HC/aeVivUR3rzLv58\nt7IafyhCXV3+FaSZDESj8PR19zJyCcqb71/L4jnnc86UCpomlWEfa7EZSZKk06QIMdqE9advYGCA\nhoYG7rvvPlasWDG0PRqNDv3+4YcfFvItj0vTBc+//kc2925i7rz+vBPwo9/8O0qT2tDjP/uHWwr2\nvk99fSVmjp7WNWDB3XcTCmVQVUE2q5BOq6TTKrmMBbfZzZrH88cGbvjWW9z0pzacNnnylyRpbI2N\njUO/+0ZeXngSCn6GcTqdNDc3s3v37kK/9GlJZTSyIovJLI6ZynnN4ilc9eofh069f/v1ldxRoGIw\nvAgI4Nt//jj1bgeZaAYdHbtqwatasTvsuLx2nr3/q+S3BnRuuM7KoWiK2ECWeCqL1azisJmpLXXh\ndcqrgyRJKoyCF4JUKsWuXbu46KKLxtynpaWl0G87poM9/ZTt7KbM5qWqKv/Kmh33XctVrz489NgC\nVFVVFeR9BywKzuzRxtYT/+s7fHSgj/6BNJoucDuseF02Sr0Odu9y8OyIqSQef2EDq978kGQuiafE\nQyqXwqyasZlsVEQUvnJBM+dOPQuLKpyEjRs3AuP7dzwZMq5TV6yxybhOzfAel9NxxoXg3nvv5fLL\nL6e2tpZDhw7xwx/+kGQyyQ033HCmL10QJxofEBz9Hl7Ia2+eePWvWLcOnn3+19Ttfh+f1cy8aZXH\n7JfK5GhtHRlRhv/a/RZ9/V04nRpmRwCbDTQNDsXgYOdByt5zG1YIJEn6dDnjQrB//36uvfZauru7\nCQaDLFq0iA0bNlBbO8YUzOMskcyQ1bJYLIJMVieb08hqOrouEEKgkX+xpvL6B0Rb67GYVaxmE6bT\nnMVTUQYHjH/xl3/L91LZY9YVPtSX4Pfv7+e26z1A/gn9pl89Rmkp9PXFURSoqgoMPdfRAb0pM27H\nWCPfkiRJp+aMC8E//uM/FiKOgktnckQTadq7YnTHEsQsKTr7kuT0HDk9hy50BHDZvfN4/X9vGjpu\nxRMv8u1f/Rlm1YxZtWA2mbCYVSwmE1aLCZvFdNLX7QsxeA/D6+9+hBCQSmcBONgbJ5Yc4ED/ATp2\n38jw0YSKpg85UkMjkfzX6++HvW0m5oamsmCGnE9IkqTC+NRdjhKJp+iKJOiND9CfjvFhuJO+ZISE\niOHLJjGbVcxWBVVVUIDktFo0NhH223jgF3+KEAK0DCktzfZ/m82Of7gL0Pn6C7djMVmxm6047Rac\ndgumkbO/HdbfPzinXFubIGz5mAHRS38yRbgnQzqjoakJTPYku34yYoEaBN965oVRX7O3F3btVJlW\n0sQFzVOOndVUkiTpNH0qCkFO0+mODtAVSRBNxulN9pHIxnE7TVjtWRx2BbPLjN97bHeKyaRw9S++\nTDAYZOSzg0VABRT+4bqnuOrXNxPPqDgyduwDdhy2wYJgswym8a8ufBgdWLr8RmJBL71anKzVSkes\nlFQa+k29DFj6wR6jrEQFbAxvDXzz8R8eE58Qg8saH2i3MivYTGvTVC77TOMx+0mSJJ2uT3Qh0DSd\njq4YXdEE0VSUvlQfupIl4LES8rgxqQq9sTQmk0L2tEaChw8jqySzaYI+F6l0kkRyAFvGhj1px2a2\nct33/xmFwYSu++XfEwem3XEl2VQF0XSEFP34qiJUBWMIIVj/zf8kf3hao3pB/i0dyaRKe7sdv91P\nS9VMLpnXyNJz6+WUEpIkFdQnthBE4yn2hqMcinfTNdCF06FSXmbD5chfh1JVFBQUxCleE/TP9/z5\niC0Kryz/Fbf8659jMavouiCVzhJJpgi3e3lqe1veO/xgQTMDcSuK0ovbnaa6JorJMnjzWsd/LQKG\nT8wn+P6ao91E2exgK2D3H71U2qtorWvhysUzmVJdHHdqS5L06fKJKwSaptPeFaOzL8qB/oMINU19\nlROrZfQZT1VFGfwGfYr3T/d+OItjLyg9mi5VVVCxEt7r4AsvvzpiYgj45fwmcnoXvvIIJbV2TOaj\n8X306x+Rf7tZkv/e2UFFwE1fj4VDnVbqfJOY5Xcwe1Ipt17xGZx2eQOZJElnxyeqEAxvBXQPdFHq\nt1LiO/4iA4MnbBVdP/Hr3/zAv9Cy7RAKsJLnRt1n5Vf+D7e8fA/xfpUPdzroTfRz/3u78k7rL4VK\nUPz7sFjDmKxeFI62Ut689WVGzifU8thy9rSb2LoNzDk/JZYyhKWE82f5mV0XkEVAkqSz6hNTCA50\n97Ovq3eoFVB3nFbAcHabilk1kzvBIMHND/wLCw4XAQFYyJDlyIIyR247UyDnJhY1sft9Oz3xCDi7\nj2kN3HZLCRZ3B1rCjKIoQ1cXpVJANkBea6DpGdr2WHCoHrwWH06nHbc7SsoKm/ZoTK+W6wxIknR2\nnd7dUuNsXzhKW7iLjyMf4/Ho1Fe5T6oIADhsZswmM9kTFIJf/OhPhn7PYEXLq5H5/Ur/+pffoDve\nh+LuIjDpABsqHOiH99rhBas3jq4P3rSWywmSSUFnJ/z3bSMHiHXKPruRxvJJzJlawdLz7XzpS7B0\nKWTUKL3JfiKJ7El9TkmSpNNV9C2CjzsjtPf0sD/WTmXQjvsUJ1tz2ExYVAu57MnXvN+zAH1oMFfQ\ncsvP2LjyDo60CpIfL8TsfQRvVRiEwpVXX048M4BN/y/M072D1xipCiJnIZG2sndAo2vdZxk5QOy7\n8SYqSjwsmG/F7z86NfWR+WAVFPTCTg4rSZJ0jKIuBHs7I+zr7uZAfwfVFQ6c9lMP12E3D3UNCcEx\nM5AOpzGYkLdYnLd97v/YzsaV+bMSRXJhfApk0xY0HYTQEaEGcgMZhKaSjPjJ9nvIanaSmh3+37ED\nxGW1XdRX2wgE8t6OeBxE1k7A4Sbozb8KSpIkqdCKtmtof1eM9p6eMyoCACZVwW4zY1LM5HLH7x76\nj881IBhZCASZjILqbedoF5FCz2O/B8Bsy2BzZLGbHWhdjSQ/nkP8w/lkOhvRY5VYdC/ZF1YycoB4\nzo+/gNNhxjvKwvb790PIHWJqpReTKu8ZkCTp7CrKQtAbS7K3q5eOWAeVQftpF4EjHDYTZtVywgHj\nf7vzYlJYRxSCDHvedzD5mz/K21c5fGJXFCidvBdfsB+z5kNNB1B0K2Z7AkfVR5RN/QiEj+GtgdCl\nf0sqBXaznYAn/xt/Tw9Euu3U+WuZUXPqC0xIkiSdqqIrBNmcxr5DEQ707ydYYj3lMYHReJwWbGYr\nqeSJP66TBCkchx8JPv8336U7kiUl+oAMR1oFApW7v9UI7e2k+90MRHxo5n7wteOevgFH439jK+2k\n/aH1jBwgrv2Tl8ikVRwWB3730UKQSsEH7yvMKJvB51umUeI+trUgSZJUaEVXCPaFo4T7uzBbNPye\nwky1XOq34TQ7SSZO5uPm73PwgEos3Yd/0n7mPZu3cABvchEbf9RFz74QiVQK3X4Ix6TtmD0RdF2h\nd8PFjBwgnvE384lEdexmB0GfC8vhG81yOdi2DSZ5pjB/SgOLmmvO7ENLkiSdpKIqBL2xJJ2RKL2p\nHirLHCc+4CSV+mw4LA6SAydzyWl+X35fKoI7FMbqSAPwI+4bevZXXM9rXEgyN4C57GPsVbtRzTkA\nNB343TPkDxAPoJoU4gkdC3ZCARcwWAS2bwe/Ws3c2mlctWSmnE9IkqRxUzRXDR3pEjoYP0B5iR2z\nuXA1qtRnw2lx0hE7/mt+sKZuxJYcJlcfzpLBhQFeu2MTATbxMlfzLi2ksXMp/0lZqBmzZ3Cfznbg\np+3ktwQABPOeXUw6bULXTFhMVqwWM+k0bN0KPrWKc6tnct3Fs7Fbi+bPIknSBFA0Z5z8LqHCXjLp\nc1tx2ezk0hZyOTCP8anf+vFfkPcN/svX48v9HkUZXASmJDv47N08ydc5sm6AQvcr34b3vwiJeobu\nQM4jcFz4BADxmAmfzYPLKQj3ZDjUbqXWPYW5tdP500vOIeApXEtIkiTpZBRF19BAKktXrL/gXUJH\nqKpCifdw91DieN1DR5+bwU70/7uazX8d5t6fvw9AjsEOnqt5CTvJw3sq8O7tkGjgyNoF+QSQY8Y1\nL5FOK2hZKz6HG7fVzc4dCuWWKZzfOJubvzBPFgFJkgxRFIWgK5KgLxnB57YUtEtouGDAjtvqIh47\nXiE4ehJfzi+Hvtt/ddMAQkDj9y7n4aYGVHKksI04buSMQ4MFwH73bOY9uxCAeMyMy+zGho+u9hJC\n9npmVTdw7edmyYnlJEkyjOGFQNN0emIDRNMRAqOsIFYotSEXPpuPSN/o/UJ79+Y/nsvmod9VoD8c\nJBGz8eTS6VQ9Ogmc+zh2bmsB6ODbSuiHMzn36QXMnDZYMFIplWzKTra3Cj1aSVNZE59pqua85hKi\niXTBPqckSdKpMnyMoCeWJJKKYbNx0hPJnY7qoBO/08tHfTYymSRWa/5J/LU7fsbw8YHP8mbe8/2H\n/KS1BPbq3aiWDKGHF9H53Q6O1tIcob+eRC6nkE1ZsVpcqIfvCtZ06O7woUTrCZU0cG71ZC6YG6Ky\nzEFXT4SuSBmhkuNPpy1JknS2GF4IBruF+igtObs3T5nNKjVBFx92e4j1xSmrGDmrZ34qOrwajbGj\nj9/6u39m0V2zMTnjQ9tCfz14rX8mp2E9cj9AVsGsmrFbDz/OmAnvrkSLVjDZN41Lzp3ChS2hoRvl\nwr0posk40bgPn1vOKyRJ0vgztGuofyBNZCBBTqTwuM5+H/mkShd+m5/oqN1D+X38X3ti3lBPvwDq\nUjksJQdHfV3rsNXHcpqCWTFjUa3EOssIfzAFra+WOudMvv75GXzpgpq8u6UDHit9qQhd0YFCfERJ\nkqRTZmiLoDeWJJKKFOwO4hOZFHLhtXvZ12VByyUxHf70bW0j9xxcW1hwtFJ22o8/cymApikoOTu5\nWBV93RXYVRfmRIC60hBfWTKZ88+tOOYYn8dKd6SfvngSTdMxmQwftpEkaYIx9KwzkM6SyiVxOsan\nHrmdFmrLvXgtfnq7j34rf++RFxk+PlB63W0AnPdAEI3By0bPe2DycV9bz9gY2F+PduAcbKl6go4K\nXCY/DZU+lsypY8m80KjHmVQFq0UhlUsxkJaL0EiSNP4MaxEIIRhIZUlraexWz7i978zJPv64P8i+\nzp5h4wT53VJ1iw9fMVRby+R7zyca0zGl2lBtA3mLx+hpJ9qAl1w8gJ5yo6fteC0BGuusuN05ejtd\nzCiv53PnVWI5zmWxdptpsBCksnhGmZZakiTpbDKsECTTOVK5NBazMnR1zXhoqPYQ8gVoj7mJRVKH\nVy4beQ/AUY5AlORAGamuerK9IRRTDhSByDhQhBmTasKmmNEwYfVqNE6OUx3ysHunh2mByVwwt5Ky\nwPEHge1WE4n+JImUbBFIkjT+DCsEiVSGVC6F3Xb2LhkdjUlVmN0YoL2vgoMd/ZhtI8cn9LxHzpII\nQleJdwXI5Rygi8GxA7OKxZbF6kxidfcQjWcpd5ZTUerlo/ed1HkbOHdqiBkNJ15TwG4z0dMnu4Yk\nSTJGwcYInn76aRoaGnA4HLS0tLB+/frj7j+QGhwfOHKZ5XiaOdlPpTeIlnSz5t5fMnx8wLfkqbx9\nFQXcwV4qmvYQavqYisZ2KqZ0EJrxEeXT2vDXdJIlicPswONwcOBjDyHHJJonVXPBvIqTmkXUbjWR\nE1mSmSyapp9wf0mSpEIqSCFYvXo1d911Fw8++CBbtmyhtbWVyy67jPb29jGPGRwoHv8WAYDFrHJe\ncxn1/jogv0Uw5dpfjXqMogrM1iwWRxqrK4XJPHhlUTqtkEpY8Vh9ZKNBfKZKZlTVsmxh5SktM2mz\nqHLAWJIkQxSkEDz55JMsX76cm2++menTp/PTn/6UyspKnnnmmTGPyeZ0snoWq8WYC5dmTvYxraac\n440PnIguINpnxmMOoEdD+E1VNFXUc2lr9SnfJW21qGS1LJmsdkrHSZIknakzPgtnMhk2bdrEsmXL\n8rYvW7aMt99+e8zjdCHQhUA1aAEWRVFYMm/kdf0n3y0jgFjEjDnnJdNTTY27ntlVU7l8yaTTWl5T\nVRSE0E+xFEmSJJ25Mx4s7u7uRtM0KiryT6rl5eV0dnaOeszGjRt5vyNKW/xjcgMWw1bjuumyr5C3\n/sCUV0gMJE58oIB4v410zI0SqWKyL0jQ4qQxlKB934enFUtPNIcl50Pr6ztmreKNGzee1muebTKu\nU1OscUHxxibjOjmNjY1ndLxhN5QNTt0gDF2Scep5Wxi8i3jwe7j/ikdIJk7wbV5Af8xOqieE0jON\nKd4GZlaGuGCOB4ft9NOpKIP5ELJJIEnSODvjFkFZWRkmk4lwOJy3PRwOU1lZOeoxLS0tqN6DqF0q\n0+s9hhWDX78ET7/wLK/s2M7s+Wm64tX0JntJxjO4PBp2hz7UXtB1SA6oxKMWcr3VuDI1zG9sYvGc\nWhbNCWI+w6khOnuSWPUAcybVU354LeMj3zpaWlrO6LULTcZ1aoo1Lije2GRcpyYajZ7R8WdcCKxW\nK/Pnz+e1117jyiuvHNr++uuvc/XVV495nKoqKIqCLsBk4DrtJV4zFS4/emqAqTUa+zrNJNJJ4pE4\nsUgGk2nwK3o2q6AMBNFjFZSYQ7RMm8oXPjuJaXXewgQiQGV8b66TJEmCAt1Qdvfdd/ONb3yDBQsW\n0NrayrPPPktnZyff+ta3xjxGVZTBAVJdgIEnv8YaOx/tr+BAJIzXn2P2VCs9USeHel0kUmmyWUjH\nHWj9blxqCY3VdcxrrGLxvAoC3sJNB6ELgaKohg2eS5I0cRWkEHz1q1+lp6eHRx55hIMHDzJ79mx+\n+9vfUltbO+YxVosJi2olndXP2vKUJ8NhU5lZ78IdaWTPvjbC++M43Q4cArS0AgNmyqx+gnVl1AXL\nWNBcxtTawndnpTMaAaflrC7OI0mSNJqCTTFx2223cdttt530/k6bBYfFQSozgGucZh8dy/RJdmZO\nr+bt95xEkwkGsgMoikLIY8NZ6qC+0ktTvY/akOuUbhI7WUII0lkdh9mB0ybXLpYkaXwZdgZ22S3Y\nzXbi6X6jQsgzvd7H1Ele+qJpuiNpTCYFr8uCz2M969NgpDM6VpMNu9UsxwgkSRp3hhUC5+FC0NWf\nMyqEY5hUhbKA/YSzhRZaKqPhMNtxOcZngR5JkqThDOuct1vNOCw2NE1B0yf2xfOptIbdbJfdQpIk\nGcKwQqAoCg6bBZvJRio9sefXSaZzg4XALguBJEnjz9ClKl12C06Lk0Ry4s64mcvpZHPIgWJJkgxj\naCEo8Tjw2/1E41nEBJ1bIdKfwWv1EvA45ECxJEmGMLZF4LDidzmxm5zEEhOvVSCEINKfIeAoIeh3\nGh2OJEkTlKGFAKDc7yJgD9AXyxgdyriLD+SwqA58DqdctF6SJMMYXggCHgc+h5dcTp1wg8Z9/WkC\n9gBBv8voUCRJmsAMLwSqqlDmc+K3++mLpY0OZ9ykMxrpNPjsHkq9DqPDkSRpAjO8EAAEDxeC+IBG\nOjMxWgVdkRQBe4AynwvTGU5hLUmSdCaK4gxks5oJ+T2UOYMc7E4aHc5ZF0tkSKdUyt1BQiVuo8OR\nJGmCK4pCAFAT9FLhKUPRbfREP71dRDlNJ9yTospTSU3QK2cblSTJcEVTCEwmlboKP1WeSnojmU9t\nF1FnTxKftYQKn18OEkuSVBSKphAAeF02KgO+T20X0ZEuoZCnnLoKv9HhSJIkAUVWCCC/i6g7kjI6\nnILJ5WSXkCRJxanoCsGRLqJqTxWRqEak/5N/o1lO09nbmaDUUS67hCRJKjrGLg02Bq/LRkOoFF0I\n9vXtxaQqeFyfzAnZNF3Q3pnAawlQ7StncmXA6JAkSZLyFGUhACgPuNB0HV3otHfvQyDwuj5ZC7fk\nNJ2O8AAO1UeNv5JptaXyngFJkopO0RYCgMpSD7ouUBSF9p596Dr4PZ+MYpDL6ezrTOA2B6gNVDKt\nphSzLAKSJBWhoi4EANVBLyaTiqoo7O3bRyqjUR6wF/WUzYlkjoNdAwTsQWr8FTTWlGAxy8FhSZKK\nU9EXAoBQiRuTqmBSTHQmwuzZH6Mq6MRpL67wdV1wqDdFPKET8tQQ8gaYWl0iu4MkSSpqxXUmPY6g\n34XLbsXVaaOrv48DhzrxuLIEi6R1cKQV4LJ4mVISoibooyLgQlGMj02SJOl4PjGFAMBptzCjrgx/\njx1nt3OodVBZ5sTlMOajaLqga1groMIboD7kx279RKVWkqQJTBHjtEZkNBodj7eRJEma0Hw+3ykf\nIzuvJUmSJjhZCCRJkia4cesakiRJkoqTbBFIkiRNcLIQSJIkTXCGFIKlS5eiqmrez3XXXWdEKDz9\n9NM0NDTgcDhoaWlh/fr1hsRxxEMPPXRMbqqqqsY9jrfeeovLL7+cmpoaVFVl1apVo8ZaXV2N0+nk\nwgsvZOfOnUUR24033nhMDltbW89qTI899hjnnXcePp+P8vJyLr/8cnbs2HHMfkbk7GRiMyJnTz31\nFHPmzMHn8+Hz+WhtbeW3v/1t3j5G5OtEcRmRq9E89thjqKrKHXfckbf9dHJmSCFQFIWbbrqJzs7O\noZ/nnntu3ONYvXo1d911Fw8++CBbtmyhtbWVyy67jPb29nGPZbimpqa83Gzbtm3cY0gkEpxzzjn8\n5Cc/weFwHHNj3OOPP86TTz7Jz372M/7whz9QXl7OJZdcQjweNzw2RVG45JJL8nI48gRTaG+++Sa3\n334777zzDm+88QZms5mLL76Yvr6+oX2MytnJxGZEzmpra3niiSfYvHkz7777LhdddBFf/vKXh/7f\njcrXieIyIlcjbdiwgZUrV3LOOefk/f+fds6EAZYuXSpuv/12I946z4IFC8Stt96at62xsVHcf//9\nBkUkxPe//30xa9Ysw95/NG63W6xatWrosa7rIhQKiUcffXRoWzKZFB6PRzz33HOGxiaEEDfccIP4\n0pe+NK5xjBSPx4XJZBL//u//LoQorpyNjE2I4siZEEKUlJSIn//850WVr+FxCWF8riKRiJgyZYpY\nu3atWLp0qbjjjjuEEGf2P2bYGMGLL75IMBhk1qxZfPe73x2Xb5LDZTIZNm3axLJly/K2L1u2jLff\nfntcYxlpz549VFdXM3nyZK699lra2toMjWektrY2wuFwXu7sdjuLFy82PHcw+I1t/fr1VFRUMH36\ndG699Va6urrGNYZYLIau6wQCg+tPFFPORsYGxudM0zRefPFFEokEra2tRZOvkXGB8bm69dZbufrq\nq1myZAli2EWfZ5IzQ+ZBuO6666ivr6eqqort27dz//33s3XrVl599dVxi6G7uxtN06ioqMjbXl5e\nTmdn57jFMdLChQtZtWoVTU1NhMNhHnnkEVpbW9mxYwclJSWGxTXckfyMlrsDBw4YEVKeSy+9lCuv\nvJKGhgba2tp48MEHueiii3j33XexWsdnGvM777yTuXPnsmjRIqC4cjYyNjAuZ9u2bWPRokWk02nc\nbje/+c1vaG5uHjpxGZWvseICY/+/Vq5cyZ49e3jhhRcA8rqFzuR/rGCF4MEHH+TRRx897j5r165l\n8eLF3HLLLUPbmpubmTJlCgsWLGDz5s3MnTu3UCF9Il166aVDv8+aNYtFixbR0NDAqlWrWLFihYGR\nnZximGTvmmuuGfq9ubmZ+fPnU1dXxyuvvMIVV1xx1t//7rvv5u2332b9+vUnlY/xzNlYsRmVs6am\nJrZu3Uo0GuWll17i+uuvZ+3atcc9ZjzyNVZczc3NhuXqgw8+4IEHHmD9+vWYTIPT2gsh8loFWhWi\nGwAAA4RJREFUYzlRzgpWCFasWMH1119/3H1qa2tH3T5v3jxMJhO7d+8et0JQVlaGyWQiHA7nbQ+H\nw1RWVo5LDCfD6XTS3NzM7t27jQ5lSCgUAgZzVVNTM7Q9HA4PPVdMKisrqampGZccrlixgn/6p39i\nzZo11NfXD20vhpyNFdtoxitnFouFyZMnAzB37lz+8Ic/8OMf/5gHHngAMC5fY8X1/PPPH7PveOXq\nnXfeobu7e6hlAoNdV+vWreO5555j+/btwOnlrGBjBKWlpUybNu24Pw6HY9Rjt23bhqZp43oCtlqt\nzJ8/n9deey1v++uvv27IpWBjSaVS7Nq1q6iKU0NDA6FQKC93qVSK9evXF1Xujujq6mL//v1nPYd3\n3nknq1ev5o033mDatGl5zxmds+PFNprxytlImqaRyWQMz9dYcY1mvHJ1xRVXsH37dt577z3ee+89\ntmzZQktLC9deey1btmyhsbHx9HN29sa2R/fRRx+Jhx9+WGzcuFG0tbWJV155RTQ1NYn58+cLXdfH\nNZbVq1cLq9Uqnn/+ebFz507xne98R3g8HrFv375xjWO4e+65R7z55ptiz549YsOGDeKLX/yi8Pl8\n4x5TPB4XmzdvFps3bxZOp1P84Ac/EJs3bx6K4/HHHxc+n0+8/PLLYtu2beKaa64R1dXVIh6PGxpb\nPB4X99xzj3jnnXdEW1ubWLNmjVi4cKGora09q7F9+9vfFl6vV7zxxhvi4MGDQz/D39OonJ0oNqNy\n9r3vfU+sW7dOtLW1ia1bt4r77rtPqKoqfve73wkhjMvX8eIyKldjWbJkSd4VmKebs3EvBO3t7WLJ\nkiWitLRU2Gw2MXXqVHHXXXeJvr6+8Q5FCCHE008/Lerr64XNZhMtLS1i3bp1hsRxxNe+9jVRVVUl\nrFarqK6uFldddZXYtWvXuMexZs0aoSiKUBRFqKo69Pvy5cuH9nnooYdEZWWlsNvtYunSpWLHjh2G\nx5ZMJsXnP/95UV5eLqxWq6irqxPLly8XHR0dZzWmkbEc+Xn44Yfz9jMiZyeKzaic3XjjjaKurk7Y\nbDZRXl4uLrnkEvHaa6/l7WNEvo4Xl1G5Gsvwy0ePOJ2cyUnnJEmSJjg515AkSdIEJwuBJEnSBCcL\ngSRJ0gQnC4EkSdIEJwuBJEnSBCcLgSRJ0gQnC4EkSdIEJwuBJEnSBCcLgSRJ0gT3/wFBopKIIzDV\nWQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x842bda0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib.patches import Ellipse\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from gif_animate import animate\n",
    "import numpy.random as random\n",
    "from DogSimulation import DogSimulation\n",
    "import filterpy.stats\n",
    "from filterpy.common import Q_discrete_white_noise\n",
    "from filterpy.kalman import KalmanFilter\n",
    "import book_format\n",
    "book_format.load_style('..')\n",
    "\n",
    "def dog_tracking_filter(R, Q=0, cov=1.):\n",
    "    dog_filter = KalmanFilter(dim_x=2, dim_z=1)\n",
    "    dog_filter.x = np.array([0, 0])   # initial state (location and velocity)\n",
    "    dog_filter.F = np.array([[1.,1],\n",
    "                             [0,1]])  # state transition matrix\n",
    "    dog_filter.H = np.array([[1.,0]])  # Measurement function\n",
    "    dog_filter.R *= R                 # measurement uncertainty\n",
    "    dog_filter.P *= cov               # covariance matrix \n",
    "    if np.isscalar(Q):\n",
    "        dog_filter.Q = Q_discrete_white_noise(2, var=Q)\n",
    "    else:\n",
    "        dog_filter.Q = Q\n",
    "    return dog_filter\n",
    "\n",
    "R = 5.\n",
    "Q = .01\n",
    "noise = 2.\n",
    "P = 20.\n",
    "dog = DogSimulation(measurement_var=R, process_var=Q)\n",
    "f = dog_tracking_filter(R=R, Q=Q, cov=P)\n",
    "random.seed(200)\n",
    "zs = []\n",
    "xs = []\n",
    "def animate_track(frame):\n",
    "    if frame > 30: return\n",
    "    \n",
    "    if frame == 0:\n",
    "        stats.plot_covariance_ellipse((0, f.x[0]), cov=f.P, axis_equal=True, \n",
    "                                      facecolor='g', edgecolor=None, alpha=0.2)\n",
    "        xs.append(f.x[0])\n",
    "        \n",
    "    z = dog.move_and_sense()[1]\n",
    "\n",
    "    zs.append(z)\n",
    "    f.update(z)\n",
    "    xs.append(f.x[0])\n",
    "    \n",
    "    stats.plot_covariance_ellipse((frame+1, f.x[0]), cov=f.P, axis_equal=True, \n",
    "                                   facecolor='g', edgecolor=None, alpha=0.5,\n",
    "                                   xlim=(-5,40), ylim=(-5,40))\n",
    "    \n",
    "    plt.plot(zs, color='r', linestyle='dashed')\n",
    "    plt.plot(xs, color='b')\n",
    "    f.predict()\n",
    "\n",
    "animate('multivariate_track1.gif', animate_track, 37, 200, figsize=(5.5, 5.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='multivariate_track1.gif'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
